{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Bayesian Analysis Using PyMC3\n",
    "\n",
    "Before we jump in to model-building and using MCMC to do wonderful things, it is useful to understand a few of the theoretical underpinnings of the Bayesian statistical paradigm. A little theory (and I do mean a *little*) goes a long way towards being able to apply the methods correctly and effectively.\n",
    "\n",
    "## What *is* Bayesian Statistical Analysis?\n",
    "\n",
    "Though many of you will have taken a statistics course or two during your undergraduate (or graduate) education, most of those who have will likely not have had a course in *Bayesian* statistics. Most introductory courses, particularly for non-statisticians, still do not cover Bayesian methods at all, except perhaps to derive Bayes' formula as a trivial rearrangement of the definition of conditional probability. Even today, Bayesian courses are typically tacked onto the curriculum, rather than being integrated into the program.\n",
    "\n",
    "In fact, Bayesian statistics is not just a particular method, or even a class of methods; it is an entirely different paradigm for doing statistical analysis.\n",
    "\n",
    "> Practical methods for making inferences from data using probability models for quantities we observe and about which we wish to learn.\n",
    "*-- Gelman et al. 2013*\n",
    "\n",
    "A Bayesian model is described by parameters, uncertainty in those parameters is described using probability distributions.\n",
    "\n",
    "All conclusions from Bayesian statistical procedures are stated in terms of *probability statements*\n",
    "\n",
    "![](images/prob_model.png)\n",
    "\n",
    "This confers several benefits to the analyst, including:\n",
    "\n",
    "- ease of interpretation, summarization of uncertainty\n",
    "- can incorporate uncertainty in parent parameters\n",
    "- easy to calculate summary statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bayesian vs Frequentist Statistics: What's the difference?\n",
    "\n",
    "Any statistical paradigm, Bayesian or otherwise, involves at least the following: \n",
    "\n",
    "1. Some **unknown quantities** about which we are interested in learning or testing. We call these *parameters*.\n",
    "2. Some **data** which have been observed, and hopefully contain information about (1).\n",
    "3. One or more **models** that relate the data to the parameters, and is the instrument that is used to learn.\n",
    "\n",
    "### The Frequentist World View\n",
    "\n",
    "![Fisher](images/fisher.png)\n",
    "\n",
    "- The data that have been observed are considered **random**, because they are realizations of random processes, and hence will vary each time one goes to observe the system.\n",
    "- Model parameters are considered **fixed**. The parameters' values are unknown, but they are fixed, and so we *condition* on them.\n",
    "\n",
    "In mathematical notation, this implies a (very) general model of the following form:\n",
    "\n",
    "<div style=\"font-size:35px\">\n",
    "\\\\[f(y | \\theta)\\\\]\n",
    "</div>\n",
    "\n",
    "Here, the model \\\\(f\\\\) accepts data values \\\\(y\\\\) as an argument, conditional on particular values of \\\\(\\theta\\\\).\n",
    "\n",
    "Frequentist inference typically involves deriving **estimators** for the unknown parameters. Estimators are formulae that return estimates for particular estimands, as a function of data. They are selected based on some chosen optimality criterion, such as *unbiasedness*, *variance minimization*, or *efficiency*.\n",
    "\n",
    "> For example, lets say that we have collected some data on the prevalence of autism spectrum disorder (ASD) in some defined population. Our sample includes \\\\(n\\\\) sampled children, \\\\(y\\\\) of them having been diagnosed with autism. A frequentist estimator of the prevalence \\\\(p\\\\) is:\n",
    "\n",
    "> <div style=\"font-size:25px\">\n",
    "> \\\\[\\hat{p} = \\frac{y}{n}\\\\]\n",
    "> </div>\n",
    "\n",
    "> Why this particular function? Because it can be shown to be unbiased and minimum-variance.\n",
    "\n",
    "It is important to note that new estimators need to be derived for every estimand that is introduced.\n",
    "\n",
    "### The Bayesian World View\n",
    "\n",
    "![Bayes](images/bayes.png)\n",
    "\n",
    "- Data are considered **fixed**. They used to be random, but once they were written into your lab notebook/spreadsheet/IPython notebook they do not change.\n",
    "- Model parameters themselves may not be random, but Bayesians use probability distribtutions to describe their uncertainty in parameter values, and are therefore treated as **random**. In some cases, it is useful to consider parameters as having been sampled from probability distributions.\n",
    "\n",
    "This implies the following form:\n",
    "\n",
    "<div style=\"font-size:35px\">\n",
    "\\\\[p(\\theta | y)\\\\]\n",
    "</div>\n",
    "\n",
    "This formulation used to be referred to as ***inverse probability***, because it infers from observations to parameters, or from effects to causes.\n",
    "\n",
    "Bayesians do not seek new estimators for every estimation problem they encounter. There is only one estimator for Bayesian inference: **Bayes' Formula**.\n",
    "\n",
    "![](images/bayes_formula.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example: Statistical hypothesis testing\n",
    "\n",
    "The *de facto* standard for statistical inference is statistical hypothesis testing. The goal of hypothesis testing is to evaluate a **null hypothesis**. There are two possible outcomes:\n",
    "\n",
    "- reject the null hypothesis\n",
    "- fail to reject the null hypothesis\n",
    "\n",
    "Rejection occurs when a chosen test statistic is higher than some pre-specified threshold valuel; non-rejection occurs otherwise.\n",
    "\n",
    "Notice that neither outcome says anything about the quantity of interest, the **research hypothesis**. \n",
    "\n",
    "Setting up a statistical test involves several subjective choices by the user that are rarely justified based on the problem or decision at hand:\n",
    "\n",
    "- statistical test to use\n",
    "- null hypothesis to test\n",
    "- significance level\n",
    "\n",
    "Choices are often based on arbitrary criteria, including \"statistical tradition\" (Johnson 1999). The resulting evidence is indirect, incomplete, and typically overstates the evidence against the null hypothesis (Goodman 1999).\n",
    "\n",
    "Most importantly to applied users, the results of statistical hypothesis tests are very easy to misinterpret. \n",
    "\n",
    "### Estimation \n",
    "\n",
    "Instead of testing, a more informative and effective approach for inference is based on **estimation** (be it frequentist or Bayesian). That is, rather than testing whether two groups are different, we instead pursue an estimate of *how different* they are, which is fundamentally more informative. \n",
    "\n",
    "Additionally, we include an estimate of **uncertainty** associated with that difference which includes uncertainty due to our lack of knowledge of the model parameters (*epistemic uncertainty*) and uncertainty due to the inherent stochasticity of the system (*aleatory uncertainty*)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Case Study: Radon contamination\n",
    "\n",
    "As a simple example, we will conduct an analysis using Gelman et al.'s (2007) radon dataset. In this dataset the amount of the radioactive gas radon has been measured among different households in all counties of several states. Radon gas is known to be the highest cause of lung cancer in non-smokers. It is believed to be more strongly present in households containing a basement and to differ in amount present among types of soil.\n",
    "\n",
    ">  the US EPA has set an action level of 4 pCi/L. At or above this level of radon, the EPA recommends you take corrective measures to reduce your exposure to radon gas.\n",
    "\n",
    "![radon](images/how_radon_enters.jpg)\n",
    "\n",
    "Let's import the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pylab as plt\n",
    "import seaborn as sns\n",
    "sns.set_context('notebook')\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>idnum</th>\n",
       "      <th>state</th>\n",
       "      <th>state2</th>\n",
       "      <th>stfips</th>\n",
       "      <th>zip</th>\n",
       "      <th>region</th>\n",
       "      <th>typebldg</th>\n",
       "      <th>floor</th>\n",
       "      <th>room</th>\n",
       "      <th>basement</th>\n",
       "      <th>...</th>\n",
       "      <th>pcterr</th>\n",
       "      <th>adjwt</th>\n",
       "      <th>dupflag</th>\n",
       "      <th>zipflag</th>\n",
       "      <th>cntyfips</th>\n",
       "      <th>county</th>\n",
       "      <th>fips</th>\n",
       "      <th>Uppm</th>\n",
       "      <th>county_code</th>\n",
       "      <th>log_radon</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5081.0</td>\n",
       "      <td>MN</td>\n",
       "      <td>MN</td>\n",
       "      <td>27.0</td>\n",
       "      <td>55735</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>N</td>\n",
       "      <td>...</td>\n",
       "      <td>9.7</td>\n",
       "      <td>1146.499190</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>AITKIN</td>\n",
       "      <td>27001.0</td>\n",
       "      <td>0.502054</td>\n",
       "      <td>0</td>\n",
       "      <td>0.832909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5082.0</td>\n",
       "      <td>MN</td>\n",
       "      <td>MN</td>\n",
       "      <td>27.0</td>\n",
       "      <td>55748</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>Y</td>\n",
       "      <td>...</td>\n",
       "      <td>14.5</td>\n",
       "      <td>471.366223</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>AITKIN</td>\n",
       "      <td>27001.0</td>\n",
       "      <td>0.502054</td>\n",
       "      <td>0</td>\n",
       "      <td>0.832909</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5083.0</td>\n",
       "      <td>MN</td>\n",
       "      <td>MN</td>\n",
       "      <td>27.0</td>\n",
       "      <td>55748</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>Y</td>\n",
       "      <td>...</td>\n",
       "      <td>9.6</td>\n",
       "      <td>433.316718</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>AITKIN</td>\n",
       "      <td>27001.0</td>\n",
       "      <td>0.502054</td>\n",
       "      <td>0</td>\n",
       "      <td>1.098612</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>5084.0</td>\n",
       "      <td>MN</td>\n",
       "      <td>MN</td>\n",
       "      <td>27.0</td>\n",
       "      <td>56469</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>Y</td>\n",
       "      <td>...</td>\n",
       "      <td>24.3</td>\n",
       "      <td>461.623670</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>AITKIN</td>\n",
       "      <td>27001.0</td>\n",
       "      <td>0.502054</td>\n",
       "      <td>0</td>\n",
       "      <td>0.095310</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5085.0</td>\n",
       "      <td>MN</td>\n",
       "      <td>MN</td>\n",
       "      <td>27.0</td>\n",
       "      <td>55011</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>Y</td>\n",
       "      <td>...</td>\n",
       "      <td>13.8</td>\n",
       "      <td>433.316718</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>ANOKA</td>\n",
       "      <td>27003.0</td>\n",
       "      <td>0.428565</td>\n",
       "      <td>1</td>\n",
       "      <td>1.163151</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 29 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "    idnum state state2  stfips    zip  region  typebldg  floor  room basement  \\\n",
       "0  5081.0    MN     MN    27.0  55735     5.0       1.0    1.0   3.0        N   \n",
       "1  5082.0    MN     MN    27.0  55748     5.0       1.0    0.0   4.0        Y   \n",
       "2  5083.0    MN     MN    27.0  55748     5.0       1.0    0.0   4.0        Y   \n",
       "3  5084.0    MN     MN    27.0  56469     5.0       1.0    0.0   4.0        Y   \n",
       "4  5085.0    MN     MN    27.0  55011     3.0       1.0    0.0   4.0        Y   \n",
       "\n",
       "   ... pcterr        adjwt  dupflag  zipflag  cntyfips  county     fips  \\\n",
       "0  ...    9.7  1146.499190      1.0      0.0       1.0  AITKIN  27001.0   \n",
       "1  ...   14.5   471.366223      0.0      0.0       1.0  AITKIN  27001.0   \n",
       "2  ...    9.6   433.316718      0.0      0.0       1.0  AITKIN  27001.0   \n",
       "3  ...   24.3   461.623670      0.0      0.0       1.0  AITKIN  27001.0   \n",
       "4  ...   13.8   433.316718      0.0      0.0       3.0   ANOKA  27003.0   \n",
       "\n",
       "       Uppm  county_code  log_radon  \n",
       "0  0.502054            0   0.832909  \n",
       "1  0.502054            0   0.832909  \n",
       "2  0.502054            0   1.098612  \n",
       "3  0.502054            0   0.095310  \n",
       "4  0.428565            1   1.163151  \n",
       "\n",
       "[5 rows x 29 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "radon = pd.read_csv('../data/radon.csv', index_col=0)\n",
    "radon.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's focus on the (log) radon levels measured in a single county (Hennepin). \n",
    "\n",
    "Suppose we are interested in:\n",
    "\n",
    "- whether the mean radon value is greater than 4 pCi/L in Hennepin county\n",
    "- the probability that any randomly-chosen household in Hennepin county has a reading of greater than 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEMCAYAAADOLq1xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxcZ33v8c9vRrusXbJW75a8L1GcOCYxScjuxAHCUgxNuJRSUqhL2nJ7uYUC5ba3LdDbFHCbXLgtDUkMlGwkOCQQEuIQ23Fiy/IuL7L2zdqtXTO/+8eMgyJka2SNfGbO/N6v13mN9MzRmd/Y0ldHz3nO84iqYowxxn08ThdgjDFmZljAG2OMS1nAG2OMS1nAG2OMS1nAG2OMS8U5XQCAiCQCVwFNgM/hcowxJlp4gUJgr6oOjX8yIgKeQLjvdLoIY4yJUhuB18Y3RkrANwHs3LmTkpISp2sxZkr+5LF9AHznY+UOV2JiTX19PRs3boRgho4XKQHvAygpKWH+/PkOl2LM1KTmBn627HvXOGjCrm27yGqMMS5lAW+MMS5lAW+MMS4VUsCLSJmI7BKRquBj6QT7zBaRn4lIpYgcE5F/FZFI6eM3xpiYE+oZ/EPANlUtA7YBD0+wz18BR1V1NbAKuBK4JyxVGmOMmbJJA15EZgPlwPZg03agXETyxu2qQJqIeIBEIAFoCGOtxhhjpiCULpQ5QIOq+gBU1ScijcH2tjH7/S/gCQLjMVOB76jqb8YfTEQygcxxzTb43RhjwiycfeQfAiqBm4A04HkR+aCq/mTcfg8AXwnj6xozJY/vqZ3S/h9dP3eGKjFmZoXSB18HFIuIFyD4WBRsH2sr8Jiq+lW1G3gGuHGC4z0ILBi3bby08o0xxlzIpAGvqq1ABbAl2LQF2K+qbeN2rQZuBxCRBOBm4NAEx+tS1TNjN6D+0t+CMcaYiYQ6iuZ+YKuIVBE4U78fQER2iMi64D4PABtF5CCBXwhVwHfDXK8xxpgQhdQHr6rHgPUTtG8a8/Ep4JbwlWaMMWY67E5WY4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxKQt4Y4xxqZACXkTKRGSXiFQFH0sn2OcREakYs/lF5O7wl2yMMSYUoZ7BPwRsU9UyYBvw8PgdVPU+VV2rqmuBjwOdwAthq9QYY8yUTBrwIjIbKAe2B5u2A+UikneRL/sk8JiqDk2/RGOMMZcilEW35wANquoDUFWfiDQG29vG7ywiCcBHgZsnOpiIZAKZ45pLplK0McaYyYUS8FP1PqBWVSsu8PwDwFdm4HWNMcaMEUrA1wHFIuINnr17gaJg+0T+APj3ixzvQeD749pKgJ0h1GKMMSZEkwa8qraKSAWwBXg0+LhfVSfqnikBNhLoornQ8bqArnFfN8WyjTHGTCbUUTT3A1tFpArYGvwcEdkhIuvG7Pdx4FlV7QhvmcYYY6YqpD54VT0GrJ+gfdO4z/8uTHUZY4yZJruT1RhjXMoC3hhjXMoC3hhjXMoC3hhjXMoC3hhjXMoC3hhjXMoC3hhjXMoC3hhjXMoC3hhjXMoC3hhjXMoC3hhjXMoC3hhjXGomFvwwxvX8fuV4Sy817X00dQ/gEeGloy2U5acxJzvF6fKMASzgjZmS+s5+frS3jqf2N1DfOfCO5z75n28C8O6yPD7xrvncsCTP1jowjrKANyYE54ZG2fbySf7fzmpG/X6uXZzL524qZVlhOl/96WH8qnzxzuW8fvIsP9hdwye+v5ebl83mnz68lozkeKfLNzHKAt6YSbxR3cHW7fto6RniniuK+fxtSyjKTH77ea9H8CJcOS+LK+dl8enrF/FnP6rg+UNN3PjNV/jY+rkUZiRf5BUuzUfXzw37MY272EVWYy5AVXn91Fk++t3dpCTE8dRn3sX/+b217wj3iSTEebh2cS6f2riQUZ+f7+2spq136DJVbcxvhRTwIlImIrtEpCr4WHqB/T4sIgdF5FDwMT+85RpzeagqPzvYxHOVTdywJI+nP3stV8zNmtIx5uWk8kfvXoRH4D93neHc0OjMFGvMBYR6Bv8QsE1Vy4BtwMPjdwiuzfpV4BZVXQlcB3SHqU5jLhtV5flDzbx+qp1rF+Xwf+9dd8n96NmpCdy3YT49AyM8uruGUZ8/zNUac2GTBryIzAbKge3Bpu1AuYjkjdv1z4BvqmozgKp2q+pgOIs15nL45dEWXjt5lmsW5rBpVSEez/RGwszJTuFD6+ZQ29HPK1VtYarSmMmFcgY/B2hQVR9A8LEx2D7WcmChiLwqIvtE5EsywRgxEckUkfljN6BkWu/CmDCprO/i5eNtrJuXxebVhWEb5riqOIMr5mTy6+NtNPfYeY+5PMJ5kTUOWA3cAlwP3AHcO8F+DwDV47adYazDmEvS0jPIk/samJudwt1ri8I+hn3TqkIS4z08ta8ev2pYj23MREIJ+DqgWES8AMHHomD7WDXAT1R1SFV7gWeAqyc43oPAgnHbxksr35jwGBrx8dieWuLjPGy5ei5xnvAPMEtNjOOu1YXUdQ6wp7oj7Mc3ZrxJv4tVtRWoALYEm7YA+1V1fGfi48CtEhAP3AQcmOB4Xap6ZuwG1E/nTRgzXS8eaaH93BAfuWrOjN6YtKYkk4W5qfzqWCvDo3bB1cysUE9T7ge2ikgVsDX4OSKyIzh6BuCHQCtwhMAvhMPA/wtvucaEX21HP7tPt7N+YQ6L8mbN6GuJCLcsz6dvaJTdp9tn9LWMCelOVlU9BqyfoH3TmI/9wJ8HN2Oigs+vPLW/nvTkeG5bfnlu25iXk0pZ/ixePdHG+gXZJMZ7L8vrmthjd7KamPbaiTZaeoa4e03RZQ3am5bm0z/s43U7izczyALexKy+oVFeqWpjWWE6ywrTL+trz8lOYWlBGjtPtDE06rusr21ihwW8iVm/rmpjeNR/2bpmxru+LI/BET8VdV2OvL5xPwt4E5O6+ofZfbqd8nlZzE5PcqSGudkpFGUmsetUO2rj4s0MsIA3MemlY60A3LR0tmM1iAjvWphLa+8Qp9r6HKvDuJcFvIk57eeG2FfTyfoF2WSmJDhay6qSDFITvOw6ddbROow7WcCbmPPaybN4PMLGsvHz5V1+8V4PVy3I5lhzLx19w06XY1zGAt7ElPZzQ7xV08kVczJJT4qMpfTWL8gBYO8Zm77AhJcFvIkpj+yqYdSvXLc41+lS3paRHE9Zfhr7azttEjITVhbwJmYMDPt4ZNcZlhakOTZy5kLK52XRMzjKqdZzTpdiXMQC3sSMn+yrp7N/hHeXOt/3Pt6ygjSS4728VdvpdCnGRSzgTUxQVR7bXcOKonTm5aQ4Xc7viPN6WDMngyONPQwM252tJjxCmmzMmGi3r7aTY829/P09q5hqN/fje2ov+nxr71BI+02mfG4Wu093UNnQ9faFV2Omw87gTUx4bHctsxLjuHtNkdOlXFBxZjL56Ynsq7FuGhMeFvDG9Tr7hnnuYBP3lBeTmhi5f7SKCGtLMqnrHKDTxsSbMLCAN673xL56hkf9fHT9XKdLmdSqkkwADjZ0O1yJcQMLeONqqsrje2q5cl4WSwsu75TAlyI7NYE5WclU1tsMk2b6Qgp4ESkTkV0iUhV8LJ1gn6+KSKuIVAS3beEv15ip2VfbxemzffzeVXOcLiVkq0oyaewepC148daYSxXqGfxDwDZVLQO2AQ9fYL9HVHVtcPtsWCo0Zhp+8lY9yfFeNq0qdLqUkK0qzkCAygY7izfTM2nAi8hsoBzYHmzaDpSLSOTdLWLMGIMjPp6rbOT2lQXMiuCLq+NlJMczLyeVyvpumyfeTEsoZ/BzgAZV9QEEHxuD7eN9REQqReRFEdkw0cFEJFNE5o/dgJJLK9+YC/vFkRZ6B0f54JXR9+21uiSDtt4hWnqsm8ZcunBeZH0IWKCqq4FvAM+IyER3azwAVI/bdoaxDmOAQPdMUUYSGxZG301DK4rSEeBwk42mMZculICvA4pFxAsQfCwKtr9NVZtVdST48S+Cz6+c4HgPAgvGbRsv9Q0YM5HWnkF2nmjj/eXFeDzidDlTlpYUz9zsFI409jhdiolikwa8qrYCFcCWYNMWYL+qto3dT0SKx3y8FpgPHJ/geF2qembsBtRf8jswZgLPVjbhV7inPPq6Z85bXpROU/egLQRiLlmoXTT3A1tFpArYGvwcEdkhIuuC+/xvETkkIgeA7wL3qmpz2Cs2JgTPHmhkeWE6i/JmOV3KJVtRlAHAkUbrpjGXJqShBap6DFg/QfumMR9/PIx1GXPJ6jr6qajr4n/cvtTpUqYlOzWBwowkDjf2cF0ETnFsIp/dyWpc57nKJgDuWh09Y98vZHlROrUd/fQOjjhdiolCFvDGdZ6rbGTNnEzmZEfevO9TtaIwAwWONvU6XYqJQhbwxlVOt53jcGMPm11w9g6Qn55IdmoCh60f3lwCC3jjKue7Z+50ScCLCCuK0jnd1mcrPZkps4A3rvJcZSNXz8+mMCPZ6VLCZkVhOj5VjrfYmHgzNRbwxjWON/dS1XKOu9a44+z9vJLsFNKS4jhsNz2ZKbKAN67xXGUjHoE7Vror4D0iLC9Mp6qllxGf3+lyTBSxgDeuoKo8e6CRDYtyyEtLdLqcsFtelM6ITznRcs7pUkwUsYA3rnC4sYcz7f1sXh25i2pPx8LcWSTFezhik4+ZKYieSbKNuYDH99Ty/KEmPALnhkZ5fE+t0yWFndcjLCtI52hTLz6/4o3CCdTM5Wdn8CbqqSoH67spnZ1GSoJ7z1mWF6UzMOLjTHuf06WYKGEBb6JeXUc/XQMjrCrJcLqUGVU6O404j3CkyUbTmNBYwJuoV9nQTZwnMNLEzRLiPJTOnsXRxh5bys+ExALeRDWfXznY0E1ZfhpJ8V6ny5lxy4vS6RoYoal70OlSTBSwgDdRbe+ZDnoHR1nt8u6Z85YUBJbys24aEwoLeBPVnj3QSLxXWFrg7u6Z82YlxjEvx5byM6GxgDdRa9Tn5/lDzSwtSCchLna+lZcXptPcM0hdR7/TpZgIF9JPhYiUicguEakKPpZeZN8lItIvIt8MX5nG/K7XT7XT0TccM90z5y0LXkx+8UiLw5WYSBfqac9DwDZVLQO2AQ9PtJOIeIPPPR2e8oy5sOcqG5mVGEdZfprTpVxWObMSKUhP4sXDtuSxubhJA15EZgPlwPZg03agXEQmWiTyC8BzQFXYKjRmAsOjfn5+qJlbl+cT742d7pnzlhWms/dMBx19w06XYiJYKD8Zc4AGVfUBBB8bg+1vE5HVwG3AP1/sYCKSKSLzx25AySXUbmLYzhNt9AyOum5q4FAtL0rHr/DSUeumMRcWllMfEYkHvgvcf/4XwUU8AFSP23aGow4TO5490EhGcjzXLZ7oD0n3K8pIoigjyfrhzUWFMnFHHVAsIl5V9QX72YuC7ecVAouAHSICkAmIiKSr6h+NO96DwPfHtZVgIW9CNDji4xdHWrhrdVFMjZ4ZS0S4ZXk+P3qzjoFhH8kJ7r/Jy0zdpD8dqtoKVABbgk1bgP2q2jZmn1pVzVXV+ao6n0CIf3eCcEdVu1T1zNgNqA/DezEx4pXjrfQN+9i8xp1TA4fq1hUFDI74efVE2+Q7m5gU6unP/cBWEakCtgY/R0R2iMi6mSrOmIk8e6CJnNQErlmY7XQpjrp6QTbpSXG8YKNpzAWENLeqqh4D1k/QvukC+391emUZM7FzQ6P88mgLH143h7gYHD0zVrzXw83L8nnpaCsjPn9MjiYyF2ffESaqvHi4maFRP+9dG9vdM+fdsaqQ7oERXj/V7nQpJgJZwJuo8kxFIyVZyVw5L8vpUiLCxtJcUhO8PH+wyelSTASygDdR4+y5IV47eZa71xQRHK0V85Livdy0LJ8XDjcz6vM7XY6JMBbwJmrsONiEz6+8d22x06VElE2rCunsH2FPdYfTpZgIYwFvosYzFY0sLUhjSUFszT0zmRuW5JGS4GWHddOYcSzgTVSo6+jnrZpO7raLq78jKd7LjUtn88LhZnx+W8rP/JYFvIkKPz3QCMDm1RbwE9m0spCz54Z5w7ppzBgW8CYq/LSikXXzspiTneJ0KRHphiV5JMV7eP6QddOY37KANxHvaFMPx1t6ee8VdnH1QlIT47ihbDbPH2rGb900JsgC3kS8ZyoaifMId66KzamBQ7VpdSFtvUO8VdvpdCkmQljAm4jm9yvPHmhkY2ku2akJTpcT0d6zdDYJcR5+VmndNCbAAt5EtLdqO2noGrCx7yGYlRjH9WV5/Ny6aUyQBbyJaE/vbyAp3sMty/OdLiUqbFpVQHPPIPvrupwuxUQAC3gTsQZHfDxX2cRtKwpITQxp4tOYd9OyfBK81k1jAizgTcR66Wgr3QMjfPBKW7I3VOlJ8Vy/JI/nKhvtpidjAW8i1xP76ilIT+Jdi3KdLiWqvHdtEa29Q+yptimEY50FvIlIbb1D/LqqjfeXF+P12MyRU3HT0nxSE7z8tKLR6VKMw0IKeBEpE5FdIlIVfCydYJ9PiEiliFSIyEER+dPwl2tixTMVDfj8ygfKrXtmqpITvNy2ooAdB5sYGvU5XY5xUKhXrh4CtqnqoyLy+8DDwHvG7fME8H1VVRFJAw6JyCuqWhnGeo3LPb6nFlXlezurKclK5o3qDptf5RLcvbaIJ/c38Ovjbdy6osDpcoxDJj2DF5HZQDmwPdi0HSgXkbyx+6lqj6qev6qTAsQDdpXHTFlT9yDNPYOUz7VVmy7VtYtzyUlN4JkD1k0Ty0LpopkDNKiqDyD42BhsfwcRuVtEDgM1wDdU9eAE+2SKyPyxG2B/h5u37avtxOsRVpdkOF1K1Ir3erhzdSG/PNJCz+CI0+UYh4T1Iquq/lRVVwBlwL0ismSC3R4AqsdtO8NZh4lePr9yoK6LZQVppCTY2PfpuKe8hKFRPztsTHzMCiXg64BiEfECBB+Lgu0TUtVa4A3grgmefhBYMG7bOLWyjVtVtfTSN+yz7pkwWFOSwaK8VJ7YV+90KcYhkwa8qrYCFcCWYNMWYL+qto3dT0SWjvk4F7gR+J0uGlXtUtUzYzfAvgMNEOiemZUYR2m+Lcs3XSLCPeUl7D3TSU17n9PlGAeE2kVzP7BVRKqArcHPEZEdIrIuuM+nReSwiFQALwHfUdUXw16xca3OvmGONfWydk6mjX0Pk3vKixGBJ/Y1OF2KcUBInZyqegxYP0H7pjEf/1kY6zIx6JmKBnyqXDE30+lSXKMwI5lrF+Xy5L56HripFI/94owpdieriQiqyg/31lGcmUxhRrLT5bjKB64spr5zgDfO2P0EscYC3kSEgw3dHGvuZd18u7gabrevKCQtMY4f773guAjjUhbwJiL8cG8dSfEe1pRY90y4JSd4uXttET872ER3v42JjyU20NjMqMf31E66z/ConyfeqmdFUTpJ8d7LUJU7hPJve15mSgJDo36eOdDAfRvmz1xRJqLYGbxx3MGGboZG/aybl+10Ka5VnJnMyuJ0tr9Rx29nFDFuZwFvHPfmmQ5yZyUyLyfF6VJc7feumsvRph4q67udLsVcJhbwxlGtPYPUdPRz1fwsRGwI30x679oikuO9/HBv6F07JrpZwBtHvVnTiUfgCpuaYMalJ8WzeU0hT+9vpHvALrbGAgt445hRv599tZ0sK0xnli2qfVnct2E+AyM+nnjLZgeJBRbwxjFHm3rpH/bZxdXLaGVxBlfMzeTR3TX4bVFu17OAN45580wHGcnxlObPcrqUmPLxDfM5fbaP106edboUM8Ms4I0jOvuHOdl6jivnZeGxi6uX1R2rCsidlcAju2qcLsXMMAt444i9wXVWr5xnF1cvt8Q4Lx+5ai4vHWuhtr3f6XLMDLKAN5fdqM/P3ppOlhakkZWS4HQ5MeneDfOI8wj//ptqp0sxM8gC3lx2hxt76BsaZf3CHKdLiVn56UlsXlPEj9+ss/lpXMwC3lx2u6vbyUlNYPFsu7jqpE9tXEj/sI/H3rC+eLeygDeXVVP3ADXt/axfkG0XVx22rDCdjaW5fP83Zxga9TldjpkBIQW8iJSJyC4RqQo+lk6wz18Hl+w7ICJvicht4S/XRLs9pzuI8wjldnE1Inxq40Jae4d4pqLR6VLMDAj1DP4hYJuqlgHbgIcn2OcN4CpVXQP8AfAjEbGleczbBkd8VNR1saYkk5QEu3M1EmwszWV5YToPvXIKn9345DqTBryIzAbKge3Bpu1AuYjkjd1PVV9Q1fNjrioBAewqmnnb/tpOhn1+rrGLqxFDRNj6nsWcPtvHc5V2Fu82oZzBzwEaVNUHEHxsDLZfyH3AKVX9nQkvRCRTROaP3YCSKVduooqqsru6g5KsZIqz7A+7SHLbigLK8mfx7V+dtOkLXCbsF1lF5HrgfwFbLrDLA0D1uG1nuOswkeX02T7aeoe4ZoGdvUcaj0fY+p5STrae4/lDzU6XY8IolICvA4pFxAsQfCwKtr+DiGwAHgXep6rHL3C8B4EF47aNUy/dRJM9p9tJjveyqiTD6VLMBDatKmRRXirf/tUJO4t3kUkDXlVbgQp+e0a+Bdivqm1j9xORq4AfAR9U1X0XOV6Xqp4ZuwE2d6mLdfUPc6Sph3Xzsoj32sjcSOT1CJ+7uYxjzb08XdHgdDkmTEL9absf2CoiVcDW4OeIyA4RWRfc51+BZOBhEakIbqvCXrGJOrtOtQOwYZF1z0Syu1YVsqo4g396sYrBERsX7wYhjVVT1WPA+gnaN435+Kow1mVcYnDExxtnOlhZnEGmzTsT0Twe4Qt3LOVj39vDD3bV8Kl3L3S6JDNNNhjZzKg3azoZGvVz3eJcp0uJeY/vCW0t1tLZs/g/v6jiw+vmkJESP8NVmZlkHaJmxoz6/Lx+6izzc1IpyUpxuhwTottXFjA44uOff1nldClmmizgzYzZcaiZrv4RNpba2Xs0KcxIZv3CbB7ZdYZDDd1Ol2OmwQLezAi/X9n2q5PkpSWypCDN6XLMFN2yrIDs1AS+9PQhGzYZxSzgzYx48UgLx1t6uXFJns0aGYWSE7x88c5lVNR18cO9v3PLi4kSFvAm7FSVb//qBPNzUlhVnOl0OeYSvW9tMRsW5vD3O45S32lL+0UjC3gTdi8fb+VwYw+fuXExXo+dvUcrEeHrH1yNX5XP/9cB66qJQhbwJqxUlX/55QmKM5N5/xXFTpdjpmlOdgpf2byC3ac7bP3WKGQBb8Jqx8FmDtR387mbSm1aApf40LoSblmez9d/fpyD9TaqJprYT6AJm+FRP19/4RhL8tP4wJU2A7RbiAj/+IHV5KUl8ukfvMnZc0NOl2RCZAFvwmb7G7XUtPfzhTuWWt+7y2SnJvDwvVfS3jfMZx/bx4jP73RJJgQW8CYsegZH+JeXTrBhYQ43LMmb/AtM1FlZnME/fmA1e6o7+NJTh1C1i66RzuaiMWHxzReO09k/zF9tWobYuHfXet8VxZxqO8e3f3WSzJR4vnDH0nf8f4c63815H10/N9wlmjEs4M20vVXTyQ921/DxDfNtQY8Y8Oe3lNE9MMLDr54mPTmez9642OmSzAVYwJtpGfH5+asnD1KQnsTnb1vidDnmMhARvrp5BT0DI3zjheMMDPv4i1vL7C+3CGQBb6blX18+xfGWXr533zpmJdq3U6zweIR/+vBakuK9fOflk3T2D/O19650uiwzjv1Emku261Q7//JSFe9dW8TNy/OdLsdcZl6P8Pf3rCIzJYGHfn2Khq4BNi7OIznB63RpJiikUTQiUiYiu0SkKvhYOsE+t4rImyIyJCLfDH+pJpK09g7ypz/cz/zcVP7u/bYyY6wSCawC9bfvW8lrJ87yb78+SWvPoNNlmaBQh0k+BGxT1TJgG/DwBPucBj4FfCNMtZkINTTq40+376d3cIR/+9iV1jVj+P1r5vH4p65hYNjHtldOsre6w4ZRRoBJA15EZgPlwPZg03agXETeMdhZVU+q6n5gdJLjZYrI/LEbYLc9RolRn5/Pba9g9+kO/v6eVTbXu3nb1Quy+ZP3lDIvO5WnKhp4dE8tfUMXjQMzw0I5g58DNKiqDyD42BhsvxQPANXjtp2XeCxzGfn9yl/+pJKfH27mr+9azvuvsN/L5p0ykuP5b9fOZ9PKAqpaevnWSyc40dLrdFkxy4k7WR8EFozbNjpQh5mCvqFRPvv4Pp7c38Bf3FLGJ69b4HRJJkJ5RLiuNI/P3LCI5AQv//H6GX56oJHhUZve4HILpfO0DigWEa+q+kTECxQF26dMVbuArrFtNn42stW29/NHP3iTqpZevnTnMgt3E5LCjGQ+e+NiXjjczK5T7VS19HJPeTELc2c5XVrMmDTgVbVVRCqALcCjwcf9qto208W5yVRv4Yap38Z9Ka9xMT6/8vqps7x0tBWvR/j4u+bzhxsXhvU1TOQKx/dTvNfDXauLWFGUwRP76vnezmquWZjD7SsKSIizqbBmWqjDH+4H/lNEvgx0AvcBiMgO4Muq+qaIXAf8EEgPPCUfAT6pqi/MQN1mBvn8yqGGbl4+3kpr7xBLC9LYvLqIrNQEp0szUWpBbip/+p5SXjwy5mzeFoSZcSEFvKoeA9ZP0L5pzMevYaNholpH3zAH6rvYW91B18AIubMSufeaeSwrTHe6NOMCCXHjzuZfq2bE5+cvb19Kqg21nRH2rxrDVJXmnkGONfdytKmH+s4BIHC2tXlNEUsK0vDY9RETZmPP5h/ZXcOvjrfy9Q+sYcOiHKdLcx0L+Bjj8yvVZ/s40tTNsaZeugZGACjJSua2FQWsKckgM8W6YszMOn82/8DNZfz3nxxgy3d38/EN8+xsPszsXzIG+FU5095HZX03hxu66Rv2Ee8VFufN4sals1lSkEZ6UrzTZZoYdPWCbH7+uXfz9ReO8f3Xz9jZfJhZwLtY+7kh9p7pYH9dF72Do8R7haUF6awqzmBJQZotim0iQnKCl69sXsEdKwvfPpu/b8M8/oedzU+b/eu5jM+vHBn80ZgAAAv5SURBVGvu4Y3qDk60nsMjsCQ/jTVzMllakG5D00zEGn82/7KdzU+bBbxLNHUP8MujLbx5poOewVHSk+K4adls1s3LJiPZul9MdLCz+fCyf7Eo5vMrr55o47HdtfzqWAuqsHj2LO5ek8OSgjS8HhsBY6LTRGfzX75rBTcvm213vk+BBXyUUVWONvXybGUjP61opKFrgNxZCXz6+kWkJsSRbTcjGZcYezb/P5+s5FOPvMmGhTl88c5lrCy2tX9DYQEf4VSV+s4BDtR38ZuT7bxa1UZD1wBej3Dt4lz+56al3Lo8cNt3uKcqMCYSXL0gm58/8G62v1HLP/+iis3feY0PlJfw+VuXUJCR5HR5Ec0C3iEjPj99Q6MMjPgC23BwC37cPTDCk/vqOX22j46+YQBmJcaxYVEOn71xMbetyCdnVuJlr9t+iZhwmsr3U5zHw5/cWMorVa08tb+Bp/c3UD4vi42Lcy/6szDVOZ3cxAJ+BrX0DHK4sZsjjT28fKyNnsERugdG6B0MBPuFCJCWFMfSwnRuWZbPypIMVhdnsLwo3YY2mpiWnODljpWFrF+QwyvHW3mrppO91R0sL0rn3aV5zMlOcbrEiGIBH0Y9gyO8WtXGr4+3sae6g9qO/refS0uKIz0pnpxZiSzITSU9OZ5ZCXEkJXhJjveSEnxMTvCSGOdBRGL6zMOYi8lOTeCe8hJuXp7PrlPt7Klu53BjD3OzU7hybhYrizNs8W8s4Ketf3iUFw4389T+RnadOsuIT8lIjufqBdnct2Eeq0syWVqYxnMHmpwu1RjXSU+K57YVBdxQlsebNZ28caaDpyoaeLaykWWF6ZTPzWLU5ycuRv/ytYC/BD6/svt0O0/sq+fnh5rpH/ZRkpXMH1y7gJuX51M+N8uGKBpzGSXGe7l2cS7vWpRDQ9cA+2o7OVDXzcGGbp6tbOTmZbO5dXkB15XmkhQfO2f2FvBTUNXSy5P7Ahd3mnsGSUuK4+41RdxTXsK6eVl4LNSNcZSIUJKVQklWCptWFnK8pZdzQ6M8f7CZH79ZT0qClxuW5HHjktlsLM1z/SgcC/hJNHcP8rODTTy1v55DDT14PcINZXl86a5l3LwsP6bOBoyJJnFeDyuKMvjo+rkMj/rZdbqdFw4384sjLew42AwEpvHYWJrLxrI81i/Idt3PswX8BE62nuOFw828eKSFA3WB5WNXl2Twlc3L2bymiFwHhicaYy5dQpyH68vyuL4sj79970qONffy6ok2dp5o45FdNXzvtWoS4jysX5DNNQtzuHJeFmtKMqP+Qm1IAS8iZcB/AjlAO3Cfqp4Yt48X+BZwO6DAP6jq98JbbvipKqfP9rGvppN9tZ3sOd3B6bN9AKyZk8l/v20Jt60oYPFsWyjYGDfweITlReksL0rn/usXMTDsY091OztPnOXVqja+8cJxAOKC+5XPzWJFUTpLCtIonZ0WVaEf6hn8Q8A2VX1URH4feBh4z7h9PgYsBkoJ/CLYLyK/VNUz4Sp2OgZHfNR3DlDb0UdNez817f1Un+2jsr6Lzv7AohfpSXGUz8viE9fO5+bl+RRmJDtctTFmpiUneLlhyWxuWDIbgM6+YfbXdfJWTWD70d66t+9bEYH5OamUzp7FnOwUijKTKc5MojgzhfyMRDKTEyJqxtZJA15EZgPlwC3Bpu3Ad0QkT1Xbxuz6e8B3VdUPtInI08CHgG+MO14mkDnuZeYB1NfXT/kN7K3u4NcnWhkZ9TPsU4ZH/Qz7/AwM+egdCtxY1DMwyrDP/46vS473UpSZzNUFs1hZnMnKogzm5aS+faF0qLOFM51TLueC2hobpvw1Z874J99pmq9hpm+gI9Cf29Zos3ZGoqn+HAEsTIKFS5L40JJCfP4CGjr7Od3Wx+mz5zjV1svhqmZe6h5iaPR3b1hMjveSnhxPelIcKQlxJMQJ8V4P8V4PCXEe4j1CQryHOPGAgCDcuDSP1SXjY3FyYzJz4j8rVPWiG3AlcHhc2xGgfFzbQeCqMZ//JfCtCY73VQJdOLbZZptttoVnu26i/HbiIuuDwPfHtSUAC4ETwIXv4XdGCbAT2AhM/U+M6BAL7xFi433GwnsEe5/neYFCYO9EXxxKwNcBxSLiVVVf8GJqUbB9rFoCXS3nX2guUDP+YKraBXRN8DpVIdRy2Y2Ze7o+Uq4nhFssvEeIjfcZC+8R7H2Oc+pCXz/p1QBVbQUqgC3Bpi3A/nH97wD/BXxKRDwikge8D3hisuMbY4yZGaFe7r0f2CoiVcDW4OeIyA4RWRfc5wfAaQLdLLuBr6nq6TDXa4wxJkQh9cGr6jFg/QTtm8Z87AP+OHylGWOMmY7IGbAZubqAv2Hi6wZuEQvvEWLjfcbCewR7nyGR4NBFY4wxLmNn8MYY41IW8MYY41IW8CESkd8XkUoRGRWRP3G6nnASkTIR2SUiVcHHUqdrCjcR+aaIVIuIishKp+uZCSKSExzZdjz4vfpkcMiyq4jI0yJyQET2i8hOEVnrdE0zSUS+cqnftxbwoasAPgI87nQhM+D8ZHJlwDYCk8m5zdPAu5ng5jsXUeDrqrpEVVcTuAHmHxyuaSZ8XFXXqOoVwDeBf3e6oJkiIuXANQRuJJ0yC/gQqeohVT0CTH3mogg2ZjK57cGm7UC52878VPU1VR1/97WrqGqHqr4ypmk3wYn83ERVu8d8moHLfibPE5FEAidcnyHwy3vKbMEPMwdoCN7HQHA6isZg+/i7lU2UEBEPgftSfup0LTNBRL4H3AoIgTUo3OhrwKOqWj1myoIpsYAPEpF9BObPmUj++QA0Jkp8GzgHfMfpQmaCqv4hgIjcS2BK8k0X/4roIiIbgKuAL0znOBbwQapa7nQNDgl1MjkTJUTkmwQW3tkcXJ/BtVT1ByLyf0UkR1Xbna4njK4HlgLnz95LgBdE5BOq+mKoB7E++Bg3hcnkTBQQkb8jsIbD+1R1yOl6wk1EZonInDGfbwY6gptrqOo/qGqRqs5X1fkEpgq+bSrhDnYna8hEZAuBPwWzgGGgD7g1eOE1qonIUgJr7mYBnQTW3D3ubFXhJSLfAu4BCoCzQLuqrnC2qvASkRXAIQJTbw8Em6tV9f3OVRVeIpIPPAOkElg7ogP4vKruc7SwGSYiZ4C7VPXQlL7OAt4YY9zJumiMMcalLOCNMcalLOCNMcalLOCNMcalLOCNMcalLOBN1BGRMyJys9N1XEik12dihwW8Mca4lAW8MRcRnLrBmKhkAW+ilogkisiDItIY3B4MTrF6/vm/FJGm4HN/GFw0YfEkx/y+iPxbcOGMPuBGEbkzuLhEj4jUichXx33NvSJSIyLtIvLFUGsUkRtEpF5E/kJEWoO1fiJ8/0Im1lnAm2j2RQKLIawF1gBXA18CEJHbgT8HbgYWE5i8KVQfBf4OSANeIzAtxX1AJnAn8Mci8r7g6ywH/g24l8AkbTkEJoaatMagAgJzmhcDnwS2iUjWFGo15oIs4E00+xjwNVVtDU6O9jcEghbgw8B/qOphVe0PPheqZ1T1N6rqV9VBVX1FVQ8GP68ksCjK+V8YHwSeU9VXg5N7/TXvXIDiYjUCjASfH1HVHQSm+F0ypX8FYy7AAt5EsyLeuQRfTbDt/HNjpzyeyvTH79hXRNaLyMsi0iYi3cD9QO5Er6OqfcDYaWsvViMEJj0bHfN5PzBrCrUac0EW8CaaNfLOJenmBtsAmnhnV8kcQjd+Br7HCayMNEdVMwisYXt+iZ2msccWkRQC3TSh1GjMjLKAN9FsO/AlEckTkVzgy8Cjwed+DHxCRJYFQ/fL03idNKBDVQdF5GoCffTn/QS4S0SuE5EEAsusjf25uliNxswoC3gTzf4WeBOoBA4C+4JtqOrzwLeAl4GTwK7g11zKIhifAb4mIr0EAvrH559Q1cPAZwmc5TcRmE+/PpQajZlpNh+8iQkisozAYhiJ4/q8jXEtO4M3riUi7xeRhOCww38EnrVwN7HEAt642aeBNuAUgeXd/hhARA6LyLkJto85Wawx4WZdNMYY41J2Bm+MMS5lAW+MMS5lAW+MMS5lAW+MMS5lAW+MMS5lAW+MMS71/wGMLVyjzxYDDwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hennepin_radon = radon.query('county==\"HENNEPIN\"').log_radon\n",
    "sns.distplot(hennepin_radon, bins=18)\n",
    "plt.axvline(np.log(4));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Bayesian Workflow\n",
    "\n",
    "The process of conducting Bayesian inference can be broken down into three general steps (Gelman *et al.* 2013):\n",
    "\n",
    "![](images/123.png)\n",
    "\n",
    "### Step 1: Specify a probability model\n",
    "\n",
    "As was noted above, Bayesian statistics involves using probability models to solve problems. So, the first task is to *completely specify* the model in terms of probability distributions. This includes everything: unknown parameters, data, covariates, missing data, predictions. All must be assigned some probability density.\n",
    "\n",
    "This step involves making choices.\n",
    "\n",
    "- what is the form of the sampling distribution of the data?\n",
    "- what form best describes our uncertainty in the unknown parameters?\n",
    "\n",
    "### Step 2: Calculate a posterior distribution\n",
    "\n",
    "The mathematical form \\\\(p(\\theta | y)\\\\) that we associated with the Bayesian approach is referred to as a **posterior distribution**.\n",
    "\n",
    "> posterior /pos·ter·i·or/ (pos-tēr´e-er) later in time; subsequent.\n",
    "\n",
    "Why posterior? Because it tells us what we know about the unknown \\\\(\\theta\\\\) *after* having observed \\\\(y\\\\).\n",
    "\n",
    "This posterior distribution is formulated as a function of the probability model that was specified in Step 1. Usually, we can write it down but we cannot calculate it analytically. In fact, the difficulty inherent in calculating the posterior distribution for most models of interest is perhaps the major contributing factor for the lack of widespread adoption of Bayesian methods for data analysis. Various strategies for doing so comprise this tutorial.\n",
    "\n",
    "**But**, once the posterior distribution is calculated, you get a lot for free:\n",
    "\n",
    "- point estimates\n",
    "- credible intervals\n",
    "- quantiles\n",
    "- predictions\n",
    "\n",
    "### Step 3: Check your model\n",
    "\n",
    "Though frequently ignored in practice, it is critical that the model and its outputs be assessed before using the outputs for inference. Models are specified based on assumptions that are largely unverifiable, so the least we can do is examine the output in detail, relative to the specified model and the data that were used to fit the model.\n",
    "\n",
    "Specifically, we must ask:\n",
    "\n",
    "- does the model fit data?\n",
    "- are the conclusions reasonable?\n",
    "- are the outputs sensitive to changes in model structure?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Probability Model\n",
    "\n",
    "Returning to our example, specifying a full probability model consists of:\n",
    "\n",
    "- a likelihood function(s) for the observations\n",
    "- priors for the unknown parameters of the model\n",
    "\n",
    "The measurements look approximately normal, so let's start by assuming a normal distribution as the sampling distribution (likelihood) for the data. \n",
    "\n",
    "$$y_i \\sim N(\\mu, \\sigma^2)$$\n",
    "\n",
    "(don't worry, we can evaluate this assumption)\n",
    "\n",
    "This implies that we have 2 unknowns in the model; the mean and standard deviation of the distribution. \n",
    "\n",
    "#### Prior choice\n",
    "\n",
    "How do we choose distributions to use as priors for these parameters? \n",
    "\n",
    "There are several considerations:\n",
    "\n",
    "- discrete vs continuous values\n",
    "- the support of the variable\n",
    "- the available prior information\n",
    "\n",
    "While there may likely be prior information about the distribution of radon values, we will assume no prior knowledge, and specify a **diffuse** prior for each parameter.\n",
    "\n",
    "Since the mean can take any real value (since it is on the log scale), we will use another normal distribution here, and specify a large variance to allow the possibility of very large or very small values:\n",
    "\n",
    "$$\\mu \\sim N(0, 10^2)$$\n",
    "\n",
    "For the standard deviation, we know that the true value must be positive (no negative variances!). I will choose a uniform prior bounded from below at zero and from above at a value that is sure to be higher than any plausible value the true standard deviation (on the log scale) could take.\n",
    "\n",
    "$$\\sigma \\sim U(0, 10)$$\n",
    "\n",
    "We can encode these in a Python model, using PyMC3 as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymc3 import Model, Normal, Uniform\n",
    "\n",
    "with Model() as radon_model:\n",
    "    \n",
    "    μ = Normal('μ', mu=0, sd=10)\n",
    "    σ = Uniform('σ', 0, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All that remains is to add the likelihood, which takes $\\mu$ and $\\sigma$ as parameters, and the log-radon values as the set of observations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "with radon_model:\n",
    "    \n",
    "    dist = Normal('dist', mu=μ, sd=σ, observed=hennepin_radon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will fit the model using a numerical approach called **Markov chain Monte Carlo (MCMC)**. This will draw samples from the posterior distribution (which cannot be calculated exactly)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [σ, μ]\n",
      "Sampling 2 chains: 100%|██████████| 4000/4000 [00:03<00:00, 1212.91draws/s]\n"
     ]
    }
   ],
   "source": [
    "from pymc3 import sample\n",
    "\n",
    "with radon_model:\n",
    "    \n",
    "    samples = sample(1000, tune=1000, cores=2, random_seed=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3wUdf7H8dfsbnolvZIQeq8BpIhUKQoIiBUL6lnRs5z6Uw/x9GyoZxcVxQYoeFKlSi/SIfQSSnoISUivuzu/Pya0E0iAJLPl83w89rFkd3bmM8Puvvf7ne/MKKqqIoQQQjgag94FCCGEEHVBAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk6IOqIoiqooSpPz/v5OUZQ39KxJCGciASeEEMIhScAJIYRwSBJwQgghHJIEnBB1y/28f/vrVoUQTkgCToi6db+iKEZFUToA/QEfRVFc9C5KCGcgASdE3fIEMoCvgYnAvUA/XSsSwkkocsFTIeqGoigq0FRV1US9axHCGUkLTgghhEOSgBNCCOGQpItSCCGEQ5IWnBBCCIckASeEEMIhmap5XvovhRBC2BKlphNKC04IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDqu5MJkI4DatV5WBmIZuP55BZUEZhmRlfdxfCfN1oGe5L2yg/PF3lIyOEvZBPq3B6OUXlzNySzPTNyWTklwHgajTg426ioKySSot2xjqTQaFrowAGtQrllk5R+Hm46Fm2EKIa1V0uR85FKRyWqqr8vDWFN38/QGG5md5NgxjRIZIejQOJ8PcAtFZddlE5e9Ly2XriNCsOnORIVhHebibu7NaQh6+PI9DbTec1EcKp1PhclBJwwikVllXy5MydrDp0im6NAnh9ZBuahfrU6LV70/L5et0xFiSk4+lq4pE+cTzYOw53F2MdVy2EQAJOiEtLyyvlge+2ciSriIk3tWJc9xgMhhp/Zs5KzCrk3SWHWLb/JLGBnrw+sg29mwbXQcVCiPNIwAlxMZn5ZYyZspH8kkq+uLszvZoGXfM81x/JZuK8vRzLLua+HrG8OKSFtOaEqDsScEL8r7ySCsZ++Sdpp0uZ+bfutIvyr7V5l1VaeHfJIb7dcJwWYT58dHtHmofVrMtTCHFFJOCEOJ/FqjLum81sO3Ga7+6Pp0eTa2+5XcyqQ1n8Y3YChWVmJt7ciru6xdTJcoRwYnLBUyHO99GKI2w8msMbt7Sps3AD6Ns8hMVPXU/3uEBenrOXf87di9lirbPlCSEuTQJOOLwNidl8svIIoztFMbZLdJ0vL9jHjW/vi+fh6+P4cVMS903bSn5pZZ0vVwhxIemiFA6tuNzMoP+sxc3FwMIJver9TCSztqbw8tw9RAd48u298cQGedXr8oVwQNJFKQTA+8sOk5ZXyjuj2+lymq2x8dH89EA3ThdXMGbKnxzKLKz3GoRwVhJwwmHtSslj2sbjjOseQ3xsgG51dIsLZPYjPTAocPtXf7I3LV+3WoRwJhJwwiGpqsqk+fsI9nbj+cHN9S6HJiHezHr4OjxdTdz59SZ2peTpXZIQDk8CTjikBbsz2JWSx3M3NsfH3TZOihwb5MUvD3fH39OVu6duZntSrt4lCeHQJOCEwymrtPDO4oO0CvdldKcovcu5QFQDT2Y9fB3BPm7cP20rh0/KPjkh6ooEnHA4P/6ZRFpeKa8Ma4nxKs4xWdfC/Nz5YXxXPFyN3PPNFtLySvUuSQiHJAEnHEpJhZkpa47Su2lQnR7Qfa2iAzz5fnxXiivMjPtmM7nFFXqXJITDkYATDuWnTUnkFFfwVP+mepdSrRZhvnxzbzxpp0sZ/91WSirMepckhEORgBMOo6TCzJdrjtG7aRBddDws4Ep0bRTAJ3d0ZHdqHv+YvZtqTrwghLgCEnDCYczckkJOcQV/H2D7rbfzDWodxotDWvD7ngy+XHtM73KEcBgScMIhmC1Wvl1/nK6xAXSOsY/W2/ke6h3HTe3CeXfJQdYePqV3OUI4BAk44RAW7c0kLa+Uh66P07uUq6IoCu+OaUezUB8mzNxJck6J3iUJYfck4ITdU1WVr9ceIy7Ii/4tQvQu56p5upr4clxnAP724zZKKyw6VySEfZOAE3Zvy/Fc9qTl80DvRhhs8Li3KxET6MXHd3Tk0MlC3vh9v97lCGHXJOCE3fthUxJ+Hi6M6mhbZy25Wn2aBfO36+OYvjmZZfsy9S5HCLslASfsWlZBGUv3ZnJr5yg8XI16l1Nrnh3YnDaRvrzw391kFZTpXY4QdkkCTti1GVuSMVtV7u4eo3cptcrVZODD2zpSWmnh2dkJWK1yfJwQV0oCTtitSouVGZuT6dMs2CGvlN0kxJuJN7Vm3ZFsvt1wXO9yhLA7EnDCbq04kEVWYbnNt95WrVpF37598fPzIzY29rLT7t+/ny5dutCgQQMaNGjAty+Pp4tfMe8uOcQRufKAEFdEAk7YrVnbUgjxcaNv82C9S7ksLy8vxo8fz+TJk6udNiIigl9//ZXc3Fyys7MZPnw4B6a/jqebkZfm7JGuSiGugAScsEvRMTHM/eELUqY+jp+vDw888AAnT55kyJAh+Pj4MGDAAE6fPg3Apk2b6NGjB/7+/rRv357Vq1efnc+0adNo2bIlPj4+xMXF8eWXX559bvXq1URFRfH+++8TEhJCeHg406ZNu+Jau3btyrhx44iLq/4gdH9/f2JjY1EUBVVVMRqNHD92lJeHtmTridP8vDXlipcvhLOSgBN2qaTcQvHBjcz/fTGHDx9mwYIFDBkyhDfffJPs7GysVisff/wxaWlpDBs2jFdeeYXc3Fzee+89Ro8ezalT2umwQkJCWLhwIQUFBUybNo2nn36aHTt2nF1OZmYm+fn5pKWl8c033/D444+fDc63334bf3//S96uhb+/P+7u7kyYMIGXXnqJMZ2juC4ukLcWH5BRlULUkAScsDuqqlJcYabD4NuJbxVHZGQkvXv3plu3bnTs2BE3NzduueUWdu7cyU8//cTQoUMZOnQoBoOBgQMH0qVLFxYtWgTAsGHDaNy4MYqi0KdPHwYNGsS6devOLsvFxYWJEyfi4uLC0KFD8fb25tChQwC8+OKL5OXlXfJ2LfLy8sjPz+fTTz+lY8eOKIrCm6PaUm62MmnBvmuatxDOQgJO2J3Nx3MxW1QGd2l+9jEPDw9CQ0Mv+LuoqIikpCRmz559Qctq/fr1ZGRkALB48WK6d+9OQEAA/v7+LFq0iOzs7LPzCQwMxGQynf3b09OToqKielhLbd/dI488wj333ENWVhaNgrx4qn9TFu3JZPn+k/VSgxD2TAJO2J1ZW1MwKNq11KoTHR3NuHHjLmhZFRcX8+KLL1JeXs7o0aN57rnnOHnyJHl5eQwdOrTG12R788038fb2vuStNlitVkpKSkhLSwPgb9fH0TzUh4nz9soFUoWohgScsCsFZZUs2puBp6sJN5fqz1xy9913s2DBApYuXYrFYqGsrIzVq1eTmppKRUUF5eXlBAcHYzKZWLx4McuWLatxLS+99BJFRUWXvJ1htVopKyujsrISVVUpKyujoqLiovNcvnw5O3fuxGKxUFBQwDPPPEODBg1o2bIlAC5GA2/c0oaM/DKmrJFrxwlxORJwwq4sSEinrNKKp1vNTssVHR3NvHnzePPNNwkODiY6OprJkydjtVrx8fHh448/ZuzYsTRo0IAZM2YwfPjwWq957dq1eHh4MHToUJKTk/Hw8GDQoEFnn2/dujXTp08HtH1vd9xxB35+fjRu3JjExESWLFmCu7v72enjYwO4uX0EX645Slpeaa3XK4SjUKrpjpGDboRNGfHpesrNVhY/1RtFse8rB1yLtLxS+r+/mgEtQ/n0zk56lyNEfarxB19acMJuJGYVkpCaz5jOUU4dbgCR/h481DuOhbsz2JuWr3c5QtgkCThhN+buTMegwPAOEXqXYhMeuj4OPw8X3l92SO9ShLBJEnDCLqiqyryENHo2CSLEx736FzgBX3cXHu4Tx6pDp9ielKt3OULYHAk4YRd2JJ8mJbeUkR0i9S7FptzXI5YgbzcmLz1U48MbhHAWEnDCLszdmY67i4Eb24TpXYpN8XQ18Xjfxmw6lsuGxBy9yxHCpkjACZtXabHy+54MBrQMxdvNVP0LnMyd3RoS4efO5GXSihPifBJwwuatO3KK3OIK6Z68BDeTkSf7NyUhJY8/DmTpXY4QNkMCTti8uTvT8fd04fpmtn3dNz2N6RxFbKAnH/5xWFpxQlSRgBM2rbjczPL9JxnWNhxXk7xdL8VkNPDoDY3Zl14g++KEqCLfGMKmLdufSWmlhZEdpXuyOiM7RhLi48aUNUf1LkUImyABJ2za3J3pRPp70LlhA71LsXluJiPjezVifWI2e1Ll7CZCSMAJm5VdVM76xGxGdIjAYHDuU3PV1J3dGuLjZmLKWmnFCSEBJ2zWwoR0LFZVuievgK+7C3d1j2HxngyScoovP/GumdpNCAclASds1txd6bQM96VZqI/epdiV8T1jMRkMfL2umuvF5adoNyEclAScsEknsovZlZLHSDmx8hUL8XVnVKdIZm9LJbuoXO9yhNCNBJywSfN2paPIlQOu2kPXx1FhsfL9xhN6lyKEbiTghM1RVZV5u9Lo1iiAcD8PvcuxS42DvenfIoSZW5IpN1v0LkcIXUjACZuzJy2fY9nFcmqua3Rvj1iyiyr4fXeG3qUIoQsJOGFz5u5Mx9VoYEjbcL1LsWu9mgTRONhLuimF05KAEzbFYlVZsDudvi2C8fNw0bscu6YoCvf2iCUhNZ+dyaf1LkeIeicBJ2zKxqPZnCosl+7JWjKqUxTebiZpxQmnJAEnbMrcnen4uJvo2yJE71IcgrebiTGdo/h9TwZZhWV6lyNEvZKAEzajrNLC0n2ZDGkThruLUe9yHMY918VQaVH5eYsc1C2ciwScsBl/HDhJUblZuidrWVywN9c3C2b65iQqLVa9yxGi3kjACZsxd2c6ob5udIsL1LsUh3NfjxhOFpSzdF+m3qUIUW8k4IRNyCupYM3hLIa3j8AoVw6odTc0CyHS30O6KYVTkYATNmHRnkwqLSojpHuyThgMCrfFR7M+MZvknBK9yxGiXkjACZswd1caTUK8aR3hq3cpDmtM5ygMCszaJq044Rwk4ITuUnJL2HI8l1s6RqIo0j1ZVyL8PejTLJjZ21Mwy2AT4QQk4ITu5iekAzC8vVw5oK7dFt+QkwXlrDl8Su9ShKhzEnBCV6qq8tuOVLrGBhAd4Kl3OQ6vf8sQgrzd+HmrdFMKxycBJ3S1L72Ao6eKGdlRBpfUBxejgTGdo1h5MIuicrPe5QhRpyTghK7m7EzD1WhgmFw5oN7cFh+NxaqyLz1f71KEqFMScEI3ZouV+QlVVw7wlCsH1JdGQV50axTA3rR8VFXVuxwh6owEnNDNxqM5nCos5xbpnqx3t3eNJq+kkpTTpXqXIkSdkYATupm7Mw1fdxM3NJcrB9S3wa3DcTUZ2J9eoHcpQtQZCTihi5IKM0v2ZTKsXbhcOUAHHq5Gmof6cORkISUVMthEOCaT3gWIepSxG/54FdJ3QmUZ+DeErg9pN4C9v8HqtyE/FVC15+MfPPf8xZjLYflEOLAQirPAIwAa94Mb/w2eAZB7HOY8Apm7Ibw93PIlNIhh2b6TPGz9hQfSDoN1HRhq4bfWH5PgwALISdT+HvE5dLzr8q/5YQRk7oWyfPDwh4bdYdAb0CBWe74kF5a+DEdXQmkueIVAy5tg4L/A5KZNU1kKq/4Ne+dA0UnwDISOd0P/f2rPz58AyZu17WpyhcguMOh1CGlZ9foyWPsu7JkNhSchIA76/h+0GlGn26RlhC970vL588+19E/9ApL+BKtZW/eb/gMx1118npumwO5fIPcYWCohqCn0eQFaDNWeT98J8yZoz8f2glumaO8F0N4L5YVw+/RrXzchqiEtOGfy813aF7VPBDQdCNmHYdFzcHyt9nx+CvhHQ/vbILY3nDp44fMXs+4D2DxF+9JqNQJQIWEGLH1Je375P7VwazYY0ndpfwPrt27jYdNCPEd8ePXhVpILqdvO/Z26DfyitZCtqYJ0aDJA+9I3mLQwmPvYueeXvqytDyq0Gqmt5+YpsO597XlVhV/uho2fgNEEHe7QQjL32Ll57PgB3Hyg7WjtPnE5/DhKCzbQttW698HgAu1vh6JMmHUvpGypQf0ZkLnn0s9fZptE+nsQ7VFKj9V3wZFlENVZW75XEBSkXXqeBxZA6WloPhRCW0HGLpg1DjIStOcXPAUFqdp7LHE5rP9Aezx5E+yfB4Pfqn69hKgF0oJzFpZK7UsHYPRU7Yvpyz7al1NesvZ4z6e02xmf94CsfXA6CRpdYr6nj2v3ncZprbbNX8Hif5yb56nD0Oh6uHUaTL8Vsg5yqrCcISn/4XDYYNo1jL+y9bBaIHEF7JoOhxZDm9EQ1UV77r6F2v0XvbTWVk08sfXcv/fP176oTyf9df16Pwfd/qYF3p+fnlu/42sg8Q8IagYPrwMX978uY/wyaNitan5J8FE7KEzXfkBEdID9c7Xnhn+stXiCW8DS/9NC785f/jo/czkc/B12zdB+sFz/Dwhre/H1u8w2MSgKN7gexqOkhMLuz+IzeGI1G6vKgEkQ2Vn7YWK1wCed4PQJOL5Oa6WfOqy9H4ZOhs+6QdZBbbrfn4Nez2g9A0LUAwk4Z2F0gW6PwqbP4LeHtG6wjAQIbQstbjo3Xep22DNL61rM2gdBzaHFsEvPt8t4OLgIdvwIRVnaF76LJ/R4Uns+uBkkroRfH4AT66HpQHb98TPxhsPkDPu+5vWfOqyFWsLPWgsnpBXc8CK0u+3qtsf5Vr+jBc7hZaAYoeeT557r9ojWClr3HqRuhcNLtdbQmW7bY2u0e1cvmNITCjMhvAMMeftc6JwJN9B+aAAoBvAJ0/5tqgrFjAQtOE7u1f7O3HthnWk7tG2w51coy9OmHfwWtL31qlc9qFI7TVru4c34JMRqrcjWI2HAa+B6iTPLRP/Pj5Iz6+Rbdaq14Gaw979QnK31EjQdCFunQmXxhdtWiDomAedMWgyDgwu0L9CTe7UvsxbDtG6zM04d1LrgQPsSbjIAXL0vPc/g5tC4LxyYrwUjaN2bZ/YvDXwdik7BoUXar/u+L9Pmi+FM97qXx0+vh0V3aNNd95i23+pipo+FI0vBJxzajtFCLbzdtW2L8+38CfKrWmRBzbU6z4jsDNHdIGn9ufVrNQIaVDVpS3K0+/Sd0HyYVuOJdVrNE7ZpwXdGeRHMfbRqfZ84F3C9n4Hfn9W6Ks907YK2P++Mz7rDqQPacrs9rG2DwMbXvOoGs3aYQEDuTtSOo1ES/4AtX2lBP+Tt6mew9CWtOzO6G7Qcrj1280faPrjDS6HJQOh4D3wzAEZ/A39+Bgkzte1yw0vQbNA1r4MQlyL74JxFSS5MH6N1rd2/BF44obUw1rwN2749N13Hu2DiaZiwQ3t+02ew8aNLz3fh01q4xT8IL2dq3Vcn1sHs+7TnAxrBA0vh5QwYv4TczT+TbXYntG0/mP8EdH0Q4sdrAzFOHb74MvKrulaDm0NoG22etenpPfBSOgx7H7IPwYyxUFGsPTf7Xi3c+r+qrV/8g9p+pIVPa897BVXV1gLumAHj5oKbr9YiTN91bhnFOfD9zZC6BTrdqw1SOSP+QXhwJfR9Bfq+DMM/uXDeZ7eBAqGttW3gF1U76+6itdI+qRzO7k6vw+A3tccPLbr866wWmPeE9mMooqPWlWqs+r0c0REeXQ8vp8Nds2D9fyCml/aDacVrMOjf2kCk2fdp+/KEqCMScM7i9AmoLNFabZGdwKOBFhigdSOBNoACtH0rgY21LyqAnKPn5nPqsHYzl2t/Zx3U7iM6gouH1uI5f57/U4P3js951Xw/A4NyQbVq00d21v6dte/itT+8Bm77CUweWihObgqz79f2wZ3pHquJ4hyt9oL0qvUt0gaJgNaiONNVW5YPxacuXL+oLtr6ndkmZ9YvtPWll3em9ZaXDN8OgvQd0OtpbV/b+ZcFMldoAzz6/AP6PK+NZgSIu+HcNM8d1kYjluXDrHu0bTDvca2L1HoNl77x1o5BNBoM/LYj9bztUdVqt1Se+z8/s5zKMvhlHOz8EeL6wr0LtffTxSRv1vYxDn7r3GCY6HitlVxZfOFgHCFqmXRROovg5tqXUOlp+H641gra86v2XMPu2v2XfaBBjDZMvCBD6xYEaNz/3Hw+q9r/8vA6rZuwYTet6+yPSdo+qmOrL5znedTFL7CEnvg27YFftIf24O/PAVVfqoFNL1670QVa3qzdirK0Ieq7ZsDM26HjOBjxqTbdug8g+4g2GhS00Ysn1kOne7Qh71u+0lqszYdpra3987Th+VFdtfA6tupcHf4xVevRTRvIMX+CFjgHf79w/VpWdVeeOggz74SKQigv0Fq/oW20ab4ZBIUZ2mjGyjJY/KL2eNtbtWDb8T3snqUN/Mk6ACmbwc1PGzxyhqunNsKx/e3aQJWEmdpt509ay6/PedOe73LbBKDhdZB7jCdc5rN0ZyZq4gEU0EaDgvZj4Mz/+QtJ2qEU85+AQ79r+w4DG8PKN7TnIztDu/P2B1otsOhZLdQbxGgDcQB+Ha/N1+gK/rEXr1uIWiAB5yxcveCuX2Hl69pghowEbaBJl/u1kYigfYEfWa59Abp4aq2V+Acv/NL6X4Pe0FqFh5dooeMRAO3vuLALDuDwUszHNzKp5F1e7xwNYeFad+bGqnAaMAnC2lS/Ht4h0GOCdkvbceFw9sQVWnfiGSmbtFtsr4sf0xXYBLyCtSCvLNPm3fFu6PPiuRbWyC9g+atacO+aoR0HF/+QVi9o3XJ3/xcWv6AFpIuntj0HvXGuy64wQ7vPT4HNX5xbflhbLeACGmk/PHbN1L70mw3W5n+pfWwNYrQBNn1e0P6vrJc5UPty2wS0HzO3fAV/vM2wgrWUqJF4DXwduj9+6XkWVK2PuUwbPHJG+zsvfK9s+1ZrJZ8Zmdt8iDbQKWEGuHjBzR+DV+CllyPENVKqOdmqnIlV1JonZuxgfWI2m1/qj5tJzl6iuzXvavd9nsdssXLd2yvpGO3PV/d00bcuIS5PqX4SjeyDE/Uiv6SSZftPMrJDpISbDTIZDYzsEMGqQ1nkFlfoXY4QtUICTtSL+QlpVJitjOlcS6P/RK0b3TmKSovKgoR0vUsRolZIwIl6MXt7Ki3DfWkT6ad3KeISWoT50ircVxtNKYQDkIATde5QZiG7U/O5VVpvNm905ygSUvM5crJQ71KEuGYScKLOzd6WgotRYaRc2NTmDW8fgdGg8NvOy5xsWQg7IQEn6lSlxcrcXWkMaBlKgJer3uWIagT7uNGnWTBzd6ZhscogamHfJOBEnVp5MIvsogpu7SLdk/ZiVKdIMvLL2HQsR+9ShLgmEnCiTs3ckkyYrzvXNw3WuxRRQwNahuLjbuK/MthE2DkJOFFn0vJKWXP4FGPjozEZ5a1mL9xdjNzULpwlezMpLr/MWVKEsHHyrSPqzC9btfMf3hYfrXMl4kqN6hRFSYWFpfsy9S5FiKsmASfqhNliZdbWFPo0CybS30PvcsQV6hLTgOgAD37bIaMphf2SgBN1YvWhU2QWlHFH14Z6lyKugqIojOoYxYaj2WTkl+pdjhBXRQJO1ImftyYT4uNGvxYhepcirtKoTpGoKszdKafuEvZJAk7Uuoz8UlYezOLWLlG4yOASuxUT6EWXmAb8tiOVaq46IoRNkm8fUetmbU3FqsLt8dI9ae9Gd47iSFYRe9MK9C5FiCsmASdqlcWq8svWZHo3DSI6wFPvcsQ1Gto2HFeTQY6JE3ZJAk7UqtWHskjPL+NOGVziEPw8XBjYKpT5CelUWqx6lyPEFZGAE7Xqx01JhPq6MaBVqN6liFoyulMkucUVrDl0Su9ShLgiEnCi1iTlFLPm8Cnu6NpQBpc4kN5NgwnyduW3ndJNKeyLfAuJWjN9czJGRZFj3xyMi9HA8PaR/LE/i/ySSr3LEaLGJOBErSirtDBrWwo3tg4j1Ndd73JELRvVKZIKi5WFe+SYOGE/JOBErViQkE5eSSV3d4/RuxRRB1pH+NI81EdO3SXsigScqBU/bUqiaYg33eMC9C5F1AFFURjVKZLtSac5nl2sdzlC1IgEnLhmCSl5JKTmM+66GBRF0bscUUdGdozEoMAcOSZO2AkJOHHNftyUhJerkVs6RupdiqhDob7u9GwSxH93pGG1yqm7hO2TgBPX5HRxBQsS0rmlUyQ+7i56lyPq2Ngu0aTllbLxaI7epQhRLQk4cU1mb0+h3GyVwSVOYmCrUPw8XJi1LUXvUoSolgScuGpWq8pPm5LpGhtAizBfvcsR9cDdxcjIDhEs2Zcpx8QJmycBJ67a2iOnSM4tYdx10npzJrd2iabCbGV+ghwyIGybBJy4aj/+mUSQtxs3tg7TuxRRj9pE+tEq3JdZ22Q0pbBtEnDiqiTlFLPyUBZ3dI3G1SRvI2cztksUe9Ly2Z8u14kTtku+mcRV+W7jCUwGRQaXOKmRHSNxNRqYvV0GmwjbJQEnrlhhWSWzt6VyU7sIOe+kk/L3dGVQ61Dm7kyj3GzRuxwhLkoCTlyx2dtSKSo3c3/PWL1LEToa2yWa0yWV/LE/S+9ShLgoCThxRSxWle82nqBLTAPaRfnrXY7QUc8mQUT4ufPz1mS9SxHioiTgxBVZeTCL5NwS7u/ZSO9ShM6MBoXb4huyPjGblNwSvcsR4i8k4MQV+Xb9cSL83LmxdajepQgbMDY+CgX4ZasMNhG2RwJO1NiBjAL+PJbDPT1iMRnlrSMg3M+Dfi1C+GVbCpUWq97lCHEB+ZYSNTZtw3E8XIzcHh+tdynChtzRtSGnCstZcUAGmwjbIgEnaiSnqJy5u9IZ3TkSf09XvcsRNqRPs2DC/dyZuUUGmwjbIgEnamTG5mQqzFbu6yGDS8SFTEYDY7tEs/bIKXxIHLEAABy3SURBVBlsImyKBJyoVoXZyg+bkujTLJgmId56lyNs0Nj4aBSQy+gImyIBJ6r1+550ThWWy4Hd4pIi/T24oXkIv2xNwSyDTYSNkIATl6WqKlPXHadxsBfXNw3Wuxxhw+7o2pCswnJWHJTBJsI2SMCJy1qfmM2+9AL+dn0cBoOidznChvVtrg02+WlTkt6lCAFIwIlqTFlzlFBfN0Z2jNS7FGHjTEYDd3VryLoj2Rw9VaR3OUJIwIlL25Oaz4bEHMb3bISbyah3OcIO3BbfEBejwo9/SitO6E8CTlzSlDVH8XE3cWe3hnqXIuxEsI8bw9qG89/t2hUnhNCTBJy4qBPZxSzem8Hd3WPwcXfRuxxhR+7pEUthuZk5O9P0LkU4OQk4cVFfrTuGyWiQQwPEFesY7U/bSD9+2HgCVVX1Lkc4MQk48RdZhWX8uj2V0Z2iCPGRK3aLK6MoCvdcF8ORrCL+PJajdznCiUnAib+YtuEElRYrf7s+Tu9ShJ26uX0EDTxdZLCJ0JUEnLhAYVklP21KYkibMBoFeeldjrBT7i5GxsZHs2z/SdLzSvUuRzgpCThxgRmbkyksM/NIn8Z6lyLs3N3dYrCqKjM2y1UGhD4k4MRZ5WYL36w/Ts8mgbSL8te7HGHnogM86d8ilJlbkik3W/QuRzghCThx1tydaWQVlkvrTdSae3vEkFNcwaI9GXqXIpyQBJwAoNJi5bNVR2kb6UevJkF6lyMcRM/GQcQFe/H9RhlsIuqfBJwAYN6udJJzS3iqf1MURU6qLGqHwaBwT/cYdqXksSslT+9yhJORgBOYLVY+XXmEVuG+9G8Zonc5wsGM7hyFj5uJqeuO6V2KcDIScIIFu9M5kVPCk9J6E3XAx92FO7s1ZPHeTFJyS/QuRzgRCTgnZ7GqfLIykRZhPgxqFap3OcJB3dczFgXtJAJC1BcJOCe3cHc6x04V82T/pnJBU1Fnwv08uLl9BL9sTSa/tFLvcoSTkIBzYtaq1luzUG8Gtw7Tu5waOXDgAP369cPPz48mTZowZ86ci0732muvoSgKf/zxx9nHZsyYQXh4OI0aNWL16tVnHz969Cg9evTAYrn0sVrfffcdvXr1+svjsbGxZ5fx3XffYTQa8fb2xtfXlw4dOrBw4UIAVq9ejcFgwNvbG29vb6Kiohg7dixbt269ms1glx7s3YjiCgszt8iB36J+SMA5sYV7MkjMKmJCP/tovZnNZkaMGMFNN91Ebm4uX331FXfffTeHDx++YLqjR4/y66+/Eh4efsFrX3zxRXbs2MEnn3zCE088cfa5J598kg8++ACj8dov6nrddddRVFREXl4eDzzwAGPHjiU3NxeAiIgIioqKKCwsZNOmTbRo0YLevXuzYsWKa16uPWgd4UfPJoFM23CcCrNV73KEE5CAc1Jmi5UPlx+meagPQ9uGV/8CG3Dw4EHS09N5+umnMRqN9OvXj549e/Ljjz9eMN0TTzzBO++8g6ur69nHcnJyiIyMJDw8nAEDBnDsmDai79dffyUyMpLu3bvXaq0Gg4Hx48dTWlp6dllnKIpCVFQU//rXv3jwwQd54YUXanXZtuzB3nGcLChn4e50vUsRTkACzkn9tiONY9nFPDOoGUY7aL0BF722mKqq7N279+zfs2fPxtXVlaFDh14wXXBwMDk5OaSmprJ8+XJat25NUVERb7zxBm+99Vat12o2m5k6dSre3t40bdr0ktONGjWKHTt2UFxcXOs12KIbmgXTNMSbr9Yek2vFiTonAeeEys0WPvzjMO2j/Oxq5GSLFi0ICQlh8uTJVFZWsmzZMtasWUNJiTb0vKioiJdeeokPP/zwL681GAx88cUXjBkzhvfee4+vv/6aiRMnMmHCBPbs2UPfvn258cYbLwjL/7Vp0yb8/f0vuCUnJ190mrCwMGbOnMmcOXPw8/O75DwjIiJQVZW8POc4CFpRFB66Po6DmYWsT8zWuxzh4Ex6FyDq34zNyaTnl/HumPZ2ddybi4sLc+fOZcKECbzzzjt06dKFsWPH4ubmBsCrr77KuHHjaNSo0UVf379/f/r37w/A7t272bZtG5MnTyY2Npb169eTkpLCgw8+yKZNmy76+u7du7N+/foLHouNja12mstJS0tDURT8/Z3n5NYjOkQweekhvl53nN5Ng/UuRzgwacE5mZIKM5+tSqR7XAA9mwTqXc4Va9euHWvWrCEnJ4elS5dy7NgxunbtCsCKFSv4+OOPCQsLIywsjJSUFMaOHcs777xzwTxUVeWJJ57g448/Jjs7G4vFQkxMDPHx8ezevbte12fOnDl06tQJLy/nufaem8nIfT1iWXv4FAczC/QuRzgwacE5mWkbTpBdVMGX45rbVevtjN27d9OsWTOsViuff/45GRkZ3HfffYAWcJWV546xio+P54MPPmDIkCEXzGPq1Kl07NiRDh06YDabKS0tZf/+/SQnJxMXV/dXMVdVlfT0dKZOncrUqVOZP39+nS/T1tzVrSGfrkxk6rrjvHdre73LEQ5KAs6J5JdU8uWao/RvEULnmAC9y7kqP/74I1OnTqWyspLevXuzfPnys12UgYEXtkiNRiMNGjTA29v77GPZ2dl89NFHbNy4EQCTycSnn35Kv379cHd3Z9q0aXVWe3p6Ot7e3qiqip+fHz169GD16tW1PoLTHvh7ujK2SxQztiTz7KBmhPt56F2ScEBKNSOZZJiTA3lj4X6+2XCc3yf0plWEr97lCL2teVe77/O8LotPyS3hhvdWc+91sUy8uZUuNQi7VOOuJ9kH5ySScor5/s8TjO0cLeEmbEJ0gCcjOkQwc0syucUVepcjHJAEnJN4e/FBXIwGnh3UTO9ShDjrsRsaU2a28N2G43qXIhyQBJwT2HI8l8V7M3mkT2NCfN31LkeIs5qE+HBjqzC+23iCwjI5CbOoXRJwDs5qVXnj9/2E+brzUO+6HyEoxJV6rG9jCsrMTN8sJ2EWtUsCzsHNS0hjd2o+zw9ujofrtZ9MWIja1i7Kn95Ng5i67jhllZe+ooMQV0oCzoGVVJh5d8kh2kb6MbJDpN7lCHFJj93QhOyicmZvS9G7FOFAJOAc2CcrE8nIL2Piza3s4nI4wnl1jwugU0N/pqw5JpfSEbVGAs5BJWYVMXXdMUZ3iiI+1j4P6hbOQ1EUnuzflLS8UmZvl1acqB0ScA5IVVVenb8XDxcj/ze0hd7lCFEjfZoF07GhP5+tTKTcLPvixLWTgHNAC3dnsCExh3/c2Jwgbze9yxGiRhRF4ekBzUjPL2PWVmnFiWsnAedgisrNvPH7ftpE+nJntxi9yxHiivRuGkSXmAZ8tuqojKgU10wCzsF8uPwwWYXlvD6ijd1cqVuIMxRF4emBzcgsKOMXacWJayQB50ASUvL4dsNxbo9vSMeGDfQuR4ir0qNxIF1jA/hsVaK04sQ1kYBzEOVmC8/NTiDU110Glgi7pigKfx/YlKzCcmbI2U3ENZCAcxAfrzjCkawi3hzVFl93F73LEeKa9GgcxHVxgXy2KpGicrPe5Qg7JQHnAPak5jNlzTFu7RxF3+YhepcjRK14YUgLcoor+GrNUb1LEXZKAs7OVZitPDc7gSBvV165SS4aKRxHh2h/hrUL5+t1x8kqKNO7HGGHJODs3Ccrj3DoZCFvjWqLn4d0TQrH8vyNzTFbrfznjyN6lyLskAScHdt6IpfPViUypnMU/VqE6l2OELUuJtCLu7rFMGtbColZRXqXI+yMBJydyi+t5O8/7yKqgSeThrfWuxwh6syEfk3wcDHy7pKDepci7IwEnB1SVZV/zt1LZkEZH93eAW83k94lCVFnAr3deKRPHMv2n2TriVy9yxF2RALODs3Zmcb8hHSeHtBUDugWTuGBXnGE+brz2oJ9WKyq3uUIOyEBZ2eSc0qYOG8fXRsF8OgNTfQuR4h64eFq5OVhLdmbVsCMLXLwt6gZCTg7UlZp4fEZO1AU+M9tHeRck8Kp3NQunB6NA5m85CA5ReV6lyPsgAScHfnXwv3sScvn/VvbE+nvoXc5QtQrRVH414jWlFRYeEcGnIgakICzE//dnsqMzck80qcxg1qH6V2OELpoEuLDA70bMWtbKtuTZMCJuDwJODtwMLOAl+fuoXtcAM8NaqZ3OULo6sl+TQn3c+efc/dhtlj1LkfYMAk4G1dQVsmjP+3A192Fj+/oiMko/2XCuXm5mfjnTa3Yn1HAdxtP6F2OsGHybWnDLFaVv/+8i+TcEj69sxMhPu56lySETRjSJowBLUN4d+khErMK9S5H2CgJOBv2zpKDrDyYxaThrenaKEDvcoSwGYqi8Oaotni5GnlmVgKV0lUpLkICzkbN3pbCV2uPMa57DOO6x+hdjhA2J8THnX/f0pbdqfl8vkouqSP+SgLOBm07kcvLc/bSs0kgE2+WS+AIcSlD24YzokMEn6w8wp7UfL3LETZGAs7GpOSW8PCP24ls4MHnd3bGRQaVCHFZ/xrehkBvV56ZtYuySove5QgbIt+eNiSnqJx7vt2C2aoy9d4u+HnK9d2EqI6fpwvvjmnPkawiJs7bi6rKuSqFRgLORpRUmBn//TbS80r55t4uNA721rskIexGn2bBPNmvCbO2pfLTZjlXpdBIwNmASouVx6fvYE9qHp/c0ZEusTJiUogr9fcBzejbPJjX5u+Ty+oIQAJOd6qq8tJve1h16BRvjGwrp+ES4ioZDAof3t6RqAYePPrTDjLzy/QuSehMAk5Hqqryr4X7mb09laf6N+XObg31LkkIu+bn4cJX93ShpMLMIz9tl0EnTk4CTieqqvLOkkNM23CC8T0b8fcBTfUuSQiH0CzUhw/GtichNY/Hpu+gwiwHgTsrCTidfLTiCFPWHOWubg35500tURS5tpsQtWVwm3DeGNmGlQezeGbWLrkKuJMy6V2AM5qy5igf/nGEMZ2jeH1EGwk3IerAXd1iKCwz8/big3i7mXhrVFv5rDkZCbh6pKoqn6xM5IPlh7m5fQTvjG6HQa7KLUSdeaRPYwrLKvls1VG83Ey8Mkx6S5yJBFw9UVWVtxYf5Ku1xxjVKZJ3R7fDKOEmRJ17blBzisstfLP+OEVlZv59Sxu57JSTkICrB1aryivz9jJjczLjusfw2vDW0nITop4oisKrN7fCx93EJysTySku55M7OuHhatS7NFHH5GdMHSs3W3h61i5mbE7m0Rsa868REm5C1DdFUXh2UHP+NaI1Kw5mcdfUTeSVVOhdlqhjEnB1KK+kgnHfbGHernSeH9ycFwa34LPPPqNLly64ublx3333XfK133//PZ07d8bX15eoqCief/55zGZz/RUvRD369NNPa/S5OF+/fv1QFOWKPhf3XBfLZ3d2Ym9aAaM+3ygXS3VwEnB1JCmnmFFfbGRXch4f3d6Bx25oAkBERASvvPIK48ePv+zrS0pK+PDDD8nOzmbz5s2sWLGC9957rz5KF6Le1fRzccb06dOv+gff0Lbh/PRgNwrKKhnx6QaW7M24qvkI2ycBVwe2J+Vyy+cbyS2u4KcHuzGiQ+TZ50aNGsXIkSMJDAy87DweffRRevfujaurK5GRkdx1111s2LChrksXQhc1/VwA5Ofn89prr/Huu+9e9fK6NgpgwYReNA314ZGfdvDOkoNyrJwDkoCrRaqq8sOfJ7j9q034uJv47dEedG1UOydOXrt2La1bt66VeQlhz1566SUeffRRwsKu7byt4X4e/PJwd+7o2pAvVh/l7qmbSc8rraUqhS2QgKslJRVmnv5lFxPn7aNXkyDmPd6TuFq65M20adPYtm0bzz33XK3MTwh7tW3bNjZs2MCECRNqZX5uJiNvjWrLu2PakZCax+AP17IgIb1W5i30JwFXC46eKuKWzzYyLyGdZwY245t74/H3dK2Vec+dO5cXX3yRxYsXExQUVCvzFMIeWa1WHnvsMT766CNMpto9wmlsl2gWPdmbuGBvJszcydO/7KKgrLJWlyHqnwTcNbBaVaZtOM6wj9eRVVjGd/d35cn+TWvtMIAlS5bw0EMPsWDBAtq2bVsr8xTCXhUUFLBt2zZuu+02wsLCiI+PByAqKop169Zd8/xjg7z49ZHr+PuApsxPSGfgB2tYvv/kNc9X6EcO9L5KqadL+Mfs3fx5LIe+zYN5e3Q7Qn3dq32d2WzGbDZjsViwWCyUlZVhMpn+8ot05cqV3HXXXcyZM4euXbvW1WoIYRNq8rnw8/MjPf1c92FKSgpdu3Zl+/btBAcH10odJqOh6sKpIbzw39089MM2hrULZ9LNrQn2cauVZYh6pKrq5W7if1SaLep3G46rrScuUVv9c7E6c3OSarVaa/z6V199VQUuuL366qtqUlKS6uXlpSYlJamqqqo33HCDajQaVS8vr7O3wYMH19VqCWe0+h3tZgNq+rk43/Hjx1VAraysrJOaKswW9ZMVh9WmLy1S201aqk7flKSaLTX/rIs6U11unb0pqnrZobEybvY8m47lMGn+Pg5mFtKrSRBvjWpLdICn3mUJcXXWVA2z7/O8vnXYuMSsIl6es4fNx3NpHeHLpOGtiY+tndHR4qrUeB+QBFwNJOeUMHnZIRYkpBPp78E/b2rJja3D5Kzkwr5JwNWYqqr8vieDf/9+gIz8Moa3j+D/hrYg3M9D79KckQRcbUjJLeHTlYn8uiMVk0HhkT6NeaRPYzlJq3AMEnBXrKTCzJTVR5my9hhGReHxvo0Z36sRnq4ynKEeScBdi6Onipi67jizt6VgMCjc2bUhj93QmJAaDCIRwm5IwF21lNwS/v37AZbsyyTI243HbmjMnd0a4u4iP37rgQTclVJVlfWJ2Xy7/jirDp3C1WjgtvhoHuvbWLohhGOSgLtmW0/k8sGyw/x5LIcwX3ce69uYWztHSy9P3ZKAq6msgjLm7Exj9vZUErOKCPJ2Y1z3GO7s1lCGBQvHJgFXazYezeaDZYfZlnQaf08X7uzakHuuiyXMT3p96oAE3OUUlFWy6mAWc3amsfbwKawqdI5pwO3x0QzvEIGbSX59CScgAVerVFVl64nTfLP+GMv2n8SoKAxuE8aYzlH0ahIkVxGvPRJw/ystr5SVB7NYti+TTcdyqLSohPu5M6pTJKM7RdXaeSOFsBsScHUmOaeE7zae4LedqeSVVBLi48bIjpEMbx9B6whfGYF9bZw74FRVJT2/jB1Jp9l4NIeNR7NJyikBIDbQkxtbhzGodSgdohtglKtrC2clAVfnys0WVh3M4tftaaw+lIXZqhLp78GAliEMbBVG10YBuJqkZXeFnCfgLFaVtNOl7EvPZ09aPnvTC9iblk9usXY5eh83E93iAujROIjeTYNoEuItv56EAAm4epZTVM6Kg1ks33+SdUdOUVZpxcPFSJfYBvRsEkSPxoG0CveVrszqOUbAWa0q+aWV5BSXk1NUQXZRBamnS0g5XUJybikpuSWkni6h0qKVaTIoNA31oW2kL20i/WgX5U+bCPt5w7z//vtMmjSJoqIivUsRTuCV67UrXryxtkLnSqrn7e3NpEmTePbZZ/UupVaUVljYkJjN+sRsNh7N5vBJ7TPv7mKgbaQfHaL9aR/tT/sof6IaeMiP8gvZVsDtTs1j6b5MrKoWWharikVVsVpVKiwqJRVmisst2n2FhZJyM3mlleQWV1z0Krv+ni40DPAkOsCThlW3VuG+NA/zsevjUCIiIsjIyNC7DOEk7CngAMLDwy842bIjySosY9OxXHYmnyYhJY+96QVUmK0ABHq50irClxZhPrQI86VluC+NQ7yceTBcjQOuXg6/P5BRwJQ12pH/igJGg4JRUTAYFFyMBrzdjHi6mvByM+Lv4UKEnzv+ni4EerkR4OVKoLcrQd7avyP8PfDzcKmPsuvds88+Ky04IS7C29vbYVpvFxPi487w9hEMbx8BQIXZyuGThexMyWNPah4HMgr54c8kyqtCz2RQaBzsTYtwH5qF+hAX5EWjYC9iA73s+kd+bbPpLkohRB2SfXB2xWyxciKnhAMZBRzMLOBgRiEHMgpIzy+7YLoIP/ezYdcoyIuoBh6E+LoT6utOiI8bLnayy+YybKuLUghhgyTgHEJRuZkT2cUcP+92LLuY46eKKCgz/2X6IG9Xgn3c8fMw4e3mgo+7CW83E97uJrxcjRgMCiaDgkFRtN42g4ICmKt2L5mtKmaL9YK/LVYVs0XFbNUet1Y9fub+/p6xdGzYoLZW2ba6KIUQNsgvWu8KRC3wdjPRJtKPNpF+FzyuqiqnSypJzyslq7CMkwXlnCzQ7rMKyigsM5OeV0pRuVm7lZmpsFivaNkGBUwGA8aqUDQaq+6rdkMZjdp9Xmllba5yjUkLTgghBACVFqs2CPC8gYBmq4qqgotRqQqyc4Fm0Oc4YumiFEII4ZBqHHB2v7dRCCGEuBgJOCGEEA5JAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDkoATQgjhkC55qq7XXnvN9NRTT9VnLUIIIcRlffTRR7FA6quvvvrXE23+j8udizLqo48+qrWihBBCiFpwHGgEnKhuwssFXGrVTOrSmULFlZHtdvVk21092XZXT7bd1bvYtkut0StVVdXtNmnSJFXP5dvrTbabbDvZdvZ1k22nz7bTe5DJazov317Jdrt6su2unmy7qyfb7upd9bar7moCQgghhF3SuwUnhBBC1AkJOCGEEA6pTgNOUZT3FEU5riiKqihKm0tMM0hRlG2KopQrivJeXdZjT2q47f6pKMo+RVESFEXZrijKjfVdpy2q4ba7X1GU3Yqi7FIUZY+iKE/Wd522pibb7bxpmyuKUiKfWU0N33OTFEXJqnrP7VIU5bP6rtMW1fR9pyjK2KrP6t6q+9Dq5l3XLbi5wPVA0mWmOQY8BEyu41rsTU223RYgXlXV9sB44BdFUTzqozgbV5Nt91+gvaqqHYAewLOKorSrj+JsWE22G4qiGIEvq6YXmhptO+AHVVU7VN0er4e67EG1205RlC7AJGCgqqptgF5AfnUzvtxxcNdMVdX1VcVdbprEqmlG1GUt9qaG227peX/uRruUeyA1PUbEQdVw2xWc96cn4AI49Yirmmy3Ki8CCwHvqpvTu4JtJ/5HDbfd08B7qqpmVr2m2nAD2QfnSO4Bjqqq6tThdiUURRmuKMo+tF+Ok1VV3aN3TbauqpV7I/AfvWuxU7dXdY0vUxTlOr2LsSOtgDhFUdYqirJDUZRXlBr8mpCAcwCKovQBXgfu0LsWe6Kq6nxVVVsDzYBxiqI017smW6YoigvwNfCIqqoWveuxQ1OARqqqtkPbJTNPUZRAnWuyFyagHTAQ6AMMAcbV5EXCjlX9CvwJGKGq6iG967FHqqomK4qyBbgJkG14aeFAY2BR1Y9nf0BRFMVXVdW/6VqZHTjTvVb17+WKoqQAbYA1+lVlN5KAX1VVLQfKFUWZB3QFfrjci6QFZ8cURYkHfgHGqKq6Q+967ImiKC3O+3cQ0BeQLsrLUFU1WVXVIFVVY1VVjQU+BL6WcKsZRVEiz/t3ByAW+UFVUzOAQYrGBegPJFT3oro+TOBjRVFSgSjgj6r9HSiKsqhqVAyKovSqmuYZ4GFFUVJluHvNth3wOeABfHne0OO2OpVsM2q47R6uOsRiF7AC+FRV1WU6lWwTarjdxEXUcNu9WTXEPQGtq3fc+a06Z1XDbfczkAXsB3YB+4Bvqp23nKpLCCGEI5IuSiGEEA5JAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDkoATQgjhkCTghBBCOKT/B5HLGe9WQ/tFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from arviz import plot_posterior\n",
    "\n",
    "plot_posterior(samples, var_names=['μ'], ref_val=np.log(4));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot shows the posterior distribution of $\\mu$, along with an estimate of the 95% posterior **credible interval**. \n",
    "\n",
    "The output\n",
    "\n",
    "    84% < 1.38629 < 16%\n",
    "    \n",
    "informs us that the probability of $\\mu$ being less than log(4) is 84% and the corresponding probability of being greater than log(4) is 16%.\n",
    "\n",
    "> The posterior probability that the mean level of household radon in Henneprin County is greater than 4 pCi/L is 0.16."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prediction\n",
    "\n",
    "What is the probability that a given household has a log-radon measurement larger than one? To answer this, we make use of the **posterior predictive distribution**.\n",
    "\n",
    "$$p(z |y) = \\int_{\\theta} p(z |\\theta) p(\\theta | y) d\\theta$$\n",
    "\n",
    "where here $z$ is the predicted value and y is the data used to fit the model.\n",
    "\n",
    "We can estimate this from the posterior samples of the parameters in the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1000/1000 [00:01<00:00, 830.53it/s]\n"
     ]
    }
   ],
   "source": [
    "from pymc3 import sample_posterior_predictive\n",
    "\n",
    "with radon_model:\n",
    "    pp_samples = sample_posterior_predictive(samples, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4633904761904762"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(pp_samples['dist'] > np.log(4)).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The posterior probability that a randomly-selected household in Henneprin County contains radon levels in excess of 4 pCi/L is ~0.46."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model checking\n",
    "\n",
    "But, ***how do we know this model is any good?***\n",
    "\n",
    "Its important to check the fit of the model, to see if its assumptions are reasonable. One way to do this is to perform **posterior predictive checks**. This involves generating simulated data using the model that you built, and comparing that data to the observed data.\n",
    "\n",
    "One can choose a particular statistic to compare, such as tail probabilities or quartiles, but here it is useful to compare them graphically.\n",
    "\n",
    "We already have these simulations from the previous exercise!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0x7b30b9d54ac8>"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAD7CAYAAAChScXIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hcV53/8feZUe+WZKtaVnHvvSRxKolJJdlNCKkbYGED7ELYZQuw/BKWh10gYfHSAwQIqZRNQthUEpPYcY1b3LuKJUuW1W11zZzfH5ITRbGtkTUzd8rn9Tz3kTRzJX1yPfPN0bmnGGstIiISHlxOBxAREd+paIuIhBEVbRGRMKKiLSISRlS0RUTCSEwgf7gxJh5YBNQCnkD+LhGRCOIe+Fhpre0b/ERAizb9BXtNgH+HiEikKgEqBj/gU9E2xlwHfAMw9HepPGCtfcaHb60FWLNmDYWFhSNKKiISraqrq1m+fPkZnxu2aBtjDPAYsNxau8sYMxtYa4x5zlrrHebbPQCFhYUUFxePLLWIiHyArzcivUD6wOcZQK0PBVtERPxs2Ja2tdYaYz4K/NEY0w6kAtcOPc8Yk0F/QR9MfSIiIn7kS/dIDPBl4CPW2rXGmAuB3xpjpltrTw069T7g/gDlFBERfOsemQvkW2vXAgx8bAemDTlvJf13OgcfZ+5JFxGR8+LL6JFqoNAYM8Vau98YMw3IBQ4PPsla2wK0DH6s/x6miIj4iy992nXGmM8AfzDGnL75+HFrbVNgo4mIyFA+jdO21j4BPBHgLCIiMgytPSLiq0sv7T9EHKSiLSISRlS0RUTCiIq2iEgYUdEWEQkjKtoiImFERVtEJIyoaIuIhBEVbRGRMKKiLSISRgK9R6RIyHpyY5XP596+pCiASUR8p5a2iEgYUdEWEQkjKtoiImFERVtEJIyoaIuIhBEVbRGRMOLLbuzFwHODHsoA0qy1mQHKJBJyntxYxRVt3QC8PsxQQQ0PlEDyZY/ICvp3ZAfAGLPSl+8TERH/G1HxNcbEAXcAKwITR0REzmWkLeYbgBpr7dahTxhjMujvOhms8HyDiYjIB420aH8C+OVZnrsPuH90cURE5Fx8LtrGmHzgEuCus5yyEvj1kMcKgTXnlUxERD5gJC3te4AXrLWNZ3rSWtsCtAx+zBhz/slEROQDRjJO+x7O3jUiIiJB4HNL21o7OZBBRERkeBpvLeJnWqdbAknT2EVEwoiKtohIGFHRFhEJIyraIiJhREVbRCSMqGiLiIQRFW0RkTCioi0iEkZUtEVEwoiKtohIGFHRFhEJIyraIiJhRAtGiQxiraXhVA/769qoauogxu0iIdZN4ZhELnDHkujpdTqiRDkVbZEBjae6eWZbDeUN7QBkJsdhraWjx8OGI428eMl9rKjdxbTuPlLi9dYRZ+iVJ1HPWsu6w428uqcOt8tw9cxcZuanMyY57t3nKxo7qH36Gf5UMJeX/3yAa2blMb8oQ7szSdCpT1uimrWWV3bX8cLOWsrGpvCFKyazfNLYdws29G+bV5KdzNd2/x+Prf8549Li+d+t1Ty5qYruPo+D6SUaqWhL1DpdsFcfbGBJSSZ3LZ1AemLsOb+ntL2BTy0v5eqZuew51sbDbx6hqb0nSIlFfCzaxpgEY8xPjDEHjTE7jTE/C3QwkUBbfeDEuwX7hjn5Pnd1uIxh+aSx/M0FxbR09vCTNw9T29oZ4LQi/XxtaX8H6AImW2tnAV8LXCSRwFt7qIFX9xxndmH6iAr2YJNzUrn3kjLcBn6xppyjTR0BSCryfsMWbWNMCnA38DVrrQWw1h4PdDCRQKlr7eILT28jOzWem+YVjOpm4rjUBD59cRmJcW4eWVv+7sgTkUDxpaVdBjQC9xtjNhtj3jDGXDT0JGNMhjGmePABFPo3rsjoeL2Wzz+9jY4eD3csLiI+xj3qn5mZHMenlpeSnhDLr9eVc+D4ST8kFTkzX4p2DFAKbLPWLgT+FXjGGJM25Lz7gPIhxxo/ZhUZtcc3VrKpvIkHbpjBuLQEv/3c9MRYPnVxKdkp8Ty2vpI9x9r89rNFBvOlaFcCfcBTANbajUADMHnIeSuBkiHHcr8lFRmlmpZOvv3SPpZPyuaWBf7/IzAlPoa/vaiU/IwEntpUxf46FW7xv2GLtrW2AfgLcCWAMWYyMA44NOS8FmttxeADqPZ/ZJGRs9by1Wd34rXwnzfNCtikmMQ4N/dcUEJuegJPbKziUP2pgPweiV6+jh65F/iKMWYn8DRwl7W2JXCxRPzrxZ11vLH/BF9aMYXxmUkB/V2JcW4+fkFxf1fJhgrdnBS/8qloW2uPWGsvtdbOstbOt9a+FOhgIv7S1evhP1/cy9TcVO65oDgovzMpPoZPXFRCRlIcj66voErDAcVPNCNSIt7PVh+hpqWT+6+fgdsVvLVCUuJj+OSFJaTGx/DrdeXUtGgCjoyeirZEtGMtnfz4jUNcMyuXZWVZQf/9aYmxfPKiEhJi3Px6XQWNp7qDnkEii4q2RLSHXtmP18KXr57mWIaMpDjuubAYay2/WlfByS6tyS3nT0VbIta+ujae3V7Dxy8sDvjNx+GMS03g7mXFnOzq5TfrK+np8zqaR8KXirZErIde2U9KfAyfvWSi01EAKMpM4mOLijjW0skftlbj7V8VQmREVLQlIm2pbOK1vfXce0kZ6UnnXm41mKblpbFiRi67alpZta/e6TgShlS0JeJYa/n2y/vJTonn4xcWOx3nA5ZPymZ+UQar9tXz6u46p+NImNF2YxLyntxY5fO5ty8p4s0DJ9hU3sQ3PjKDpLjQe4kbY7hxbgHH27r50u/f4YW8NMf73CV8qKUtEcXrtXzn5f2Mz0zk1kVFTsc5qxi3i9sWF2GBv39yq25Mis9CrxkiMgov7KxlT20bK2+dS1xMaLdJMpPjePDm2dz7+FYefGUfX712+jnPH+lfHBKZQvtVLTICHq/lu6/uZ2puKjfMyXc6jk8+PDOPO5cW8Yu3ytlU3uR0HAkDKtoSMbZUNlPR2ME/r5iCK4jT1Ufry1dPo3BMIv/8h3fo6OlzOo6EOBVtiQi9Hi+r9h1nwYQxXD51nNNxRiQ5PoYHb55DZWMH33ppn9NxJMSpaEtEWH+4kbauPv71w1MDtlZ2IC0tzeKeC4r5zfpKtlY1Ox1HQpiKtoS9zh4Pbx44weScFBaXZDod57x9acUUctLi+dpzu/B4NVtSzkxFW8LemkMn6Oz1cNX0XKejjEpKfAxfu246u4+18cTGSqfjSIhS0Zaw1tbVy9pDDcwuTCc/I9HpOKN27aw8LpqYzYOv7OfESS3jKh+koi1h7fW99Xi9cOW0HKej+IUxhq9/ZAZdvR6+99oBp+NICPKpaBtjKowx+4wx2weOFYEOJjKc+pNdbKlsYnFJJlkp8U7H8ZuysSncsWQCv337qDYGlg8YSUv7Zmvt3IHjlYAlEvHRq7uPE+t2cVmYDfHzxT9cPpHEWDffeVlDAOX9/DaN3RiTAWQMebjQXz9fIstIpmSfSWVjO3tq2/jQtBxS4iNvNYaslHg+c2kZD76yn7crmlhUHL6jYsS/RtLSfsIYs8MY8+OBAj3UfUD5kGONHzKKvI+1lpd21ZGaEMNFE7OdjhMwn7iwhJy0eL710j6sNkyQAb4W7eXW2jnAIsAAPzzDOSuBkiHHcn+EFBlsb20bVU0dXDE1J+QXhRqNxDg3/3D5JLZUNvPWoQan40iI8OkVb609OvCxG/gxcOEZzmmx1lYMPoBqf4YV8Xgtr+w+ztiUeBZMGON0nIC7ZWEh+ekJrHztoFrbAvhQtI0xycaY9IHPDfAxYHugg4mcyZbKZk6c6mbFjBzcYbQo1PmKj3Hz2csmsqWymUMnNJJEfGtp5wBvGGN2ALuAycBnA5pK5Aw6ezy8uqeOCVlJTMtLczpO0Jxubb++t16tbRm+aFtrj1hr51lrZ1trZ1hrb7HW1gYjnMhgr+09TmePh+tn54flolDnKz7GzWcum0hVUweHT7Q7HUccFrl3cSSi1LZ2suFII0tKMyNiuvpI3bKgkJT4GNYcPOF0FHGYiraEPGstf3rnGIlxbj4UIdPVRyoh1s2FZVkcrD/FsZZOp+OIg1S0JeS9U91KRWMHK2bkhuTu6sGyuCSL+BgXq9Xajmoq2hLSuns9vLSrlsIxiVExxO9cEuPcLC7JZGd1K03tPU7HEYeoaEtIW7W/npNdfVw/Ox9XFN18PJsLy7JxuYwm20QxFW0JWfUnu1h7qIEFE8YwPjPJ6TghIS0xltkF6Wytaqar1+N0HHGAiraEJK+1PLethrgYFytmhPeONP52QVk2PX1etlRqL8lopKItIentiiYqGju4dlZeRK7iNxoFYxIpykxi/ZFGvJpsE3VUtCXktHb28vKuOsrGJjO/KLpvPp7NBWVZNLX3cKDupNNRJMhUtCWkWGv54/YavNZy07zCqJr5OBIz8tNJS4hh3ZFGp6NIkKloS0jZWdPKvrqTfGhaDpnJcU7HCVlul2FpaRaH6k9xvK3L6TgSRCraEjI6uvv4045aCjISuaAscjc38JeFxZnEuAwb1NqOKiraEjJe3FVLZ08ffzW/ICqWXR2tlPgY5hRmsLWqmc4eDf+LFiraEhIOHj/J1qoWLp40lrz06FsQ6nwtK8ui12PZXNnkdBQJEhVtcVx3n4fntteQnRIfkTurB1J+RiLFWUls0PC/qKGiLY57bc9xmjt6uWleAbFuvSRHallZNs0dveyr1fC/aKB3iDjqaFMH6w43sqQkk5LsZKfjhKXpeWmkJcSwoVw3JKPBiIq2MeZ+Y4w1xswMVCCJHh6v5dltNaQmxGiq+ii4XYbFJZkcqj9Fw8lup+NIgPlctI0x84GlQFXg4kg0WXe4gbq2Lq6fk09CrNvpOGFtUXEmLgMb1dqOeD4VbWNMPPAj+jf01d0OGbXmjh5e23ucqbmpTI+iTXoDJTUhlhn56Wypaqanz+t0HAkgX1fi+Q/gcWtt+dmmFRtjMoCMIQ8XjiKbRKjT24cBXD8nujbpDaSlpVnsrGnlneoWp6NIAA1btI0xy4BFwL8Nc+p9wP3+CCWR7ZXddeyrO8nVM3MZk6Sp6v5SnJVETlo8G440Yq3V/wwjlC/dI5cAU4FyY0wF/a3nV4wxVw05byVQMuRY7r+oEglOdffxwPN7yE1L0FR1PzOmfz2S2tYutlaptR2phi3a1tpvWWvzrbXF1tpioBpYYa19dch5LdbaisHHwLki7/ruq/s5frKLG+dpqnogzB2fQXyMi8c3VDodRQJE47QlaHbVtPLougpuX1xEkbYPC4j4GDfzisbwwo5aGk5p+F8kGnHRHmhx7wpEGIlc1lq+/qfdjEmK418+PNXpOBFtaUkmPR4vv9t81OkoEgBqaUtQvLSrjrcrmvnHqyaTnhjrdJyINi4tgWWlWTyxoQqPVyN0I42KtgRcV6+H/3ppL1NyUrl14Xin40SFu5dNoKalk1X76p2OIn6moi0B9+t1FRxt6uTfr5tGjBaECooPTc8hJy2ex3RDMuLoHSQB1drRy4/+cojLpoxl+aSxTseJGrFuF7cvnsDqAycob2h3Oo74kYq2BNTDqw9zsqtPNx8dcNvi8cS4DE+otR1RVLQlYE6c7OZXayu4YU4+07S+SNCNS0tgxcxcfrf5qLYjiyC+rj0iUerJjb4v6nj7kqL3ff2jvxyix+Pli1dO9ncsGcbpf7e89ATauvr46rM7WVicecZzh/67SWhTS1sCora1kyc3VnHLgkJtbuCgkqxkxqXGs6G8fz0SCX9qaUtAPPzmEbzW8veXTwzq7x3JXwbR4PR6JM+/c4zq5k7GayZq2FNLW/yu4VQ3T79dxU3zCigcoyLhtHkD65FsOKINEiKBirb43a/WltPd5+XeS8ucjiJAfKybeUUZ7Khppb27z+k4Mkoq2uJXbV29/GZdJVfPzKVsbIrTcWTAkpIsPF7L5spmp6PIKKloi189tr6Sk919fPbS4PZly7nlpCVQkp3MpvJGvLohGdZUtMVvevq8PLquguWTsplZkO50HBliaWkWzR29HKg76XQUGQUVbfGbl3bVUn+ym09cVOJ0FDmD6XlppCbEsEE7toc1FW3xC2stj7xVTunYZC7RGiMhye0yLCrO5MDxUzRqg4SwpXHa4hdVTR3sqG7lGzfOxKVtxHwW7HHli4szeWN/PRvLm7hmVl5Qf7f4h1ra4hdrDzeSlhDDX88vcDqKnENaYizT89PZUtlMr8frdBw5Dz4VbWPMc8aYd4wx24wxa4wxcwMdTMJHa2cve461ctviIpLi9MdbqFtakklnr4cd1dqxPRz52tL+G2vtHGvtPOAh4JcBzCRhZnNlE14LdyyZ4HQU8UFJ9sB6JEeanI4i58Gnom2tbR30ZTqgv6sEoH/CRkUzk8alUJSlKevhwBjDktIsalo6OdrU4XQcGSGf/5Y1xvwCuAowwIfP8HwGkDHk4cJRpZOQd+D4SVo7e7l+tm5qhZN54zN4ZXcdGzX8L+z4XLSttX8LYIy5C3gQuGbIKfcB9/svmoSDjeX9NyCn5GqTg3CSEOtm3vgMtlQ209TeQ2ZynNORxEcjHj1irX0MuMwYkzXkqZVAyZBj+agTSshqau/h4PFTLCzOxK1hfmFnSWkWfV7L7zcfdTqKjMCwRdsYk2KMGT/o6+uBpoHjXdbaFmttxeADqPZ3YAkdWwYWH1o4YYzDSeR85KYlUJyVzGMbKvF4tR5JuPClpZ0M/N4Ys9MYsx34InC91TYYUc1rLduqmpk4LoWMJP1pHa6WlWVR3dzJqn31TkcRHw3bp22tPQ4sDUIWCSPlDe20dPayYkau01FkFKbnpZGXnsCj6yq4cnqO03HEB5oRKedla2Uz8TEupufrBmQ4c7sMdy6dwFuHGjhUr9X/woGKtoxYd6+HXcdamV2YQaxbL6Fwd+ui8cS5XTy6rtLpKOIDzTmWEdt1rJVej2VB0dBh+RKOXt19nBn5afz27aOUZCeTEOs+67m3LykKYjI5EzWTZMS2VLaQnRKnnb0jyLKyLHo83ndHBEnoUtGWEWlq76GisZ35RWMwRmOzI0XhmCTGj0lkwxFtRxbqVLRlRLZWNWOAuePVNRJplpVl09jew6H6U05HkXNQ0RafnR6bXaax2RFpZkEaqfExrD+s9UhCmYq2+KyioZ3mjl7mF2kGZCSKcblYVJLJgeMntR1ZCFPRFp9trRoYm52nsdmRanFJJsbAhiNqbYcqFW3xSXefh101bcwqSCcuRi+bSJWWEMvMgnQ2VzbT3edxOo6cgd594pPdNW30eLws0OJQEW9ZaRbdfV62VWk7slCkoi0+2VLVTFZyHEUamx3xijKTyM9IYMORRrQuXOhR0ZZhNbX3UN7QzvwJGpsdDYwxLCvNpv5kN4dPtDsdR4ZQ0ZZhbTvaPzZ7nsZmR43ZhekkxblZrxuSIUdFW87Jay1bK5spHZussdlRJNbtYlFxJvtq22hu73E6jgyioi3nVNnYobHZUWrJwPA/bf4bWlS05Zy2VjYTF+NiRn6601EkyDKS4piWl8bbFc30erxOx5EBvuwRmWWMedEYs98Ys8MY84wxZmwwwomz2rv72FnTymyNzY5ay0qz6Oz18M5RDf8LFb68Ey3wHWvtFGvtbOAw8K3AxpJQ8MLOWo3NjnIl2cnkpMWzXsP/QsawRdta22StfWPQQxuACQFLJCHj95uPkp0Sr7HZUez08L/a1i4qGjucjiOMsE/bGOMCPgM8f4bnMowxxYMPoNAvKSXojpw4xdsVzSzU2OyoN3d8BgmxLg3/CxEj7aj8AXAK+OEZnrsPKB9yrBlVOnHMH7ZU43YZ5mpLsagXF+Ni4YRM9hxrpba10+k4Uc/nom2MeQiYBNxqrT3TreSVQMmQY7k/Qkpw9Xm8/O/Wai6dPJa0hFin40gIWFqahbXwxIYqp6NEPZ+KtjHmm8AC4EZr7RkX2rXWtlhrKwYfQLX/okqwrDnYwPG2bm5ZON7pKBIiMpPjmJKbylObqujq1ep/TvJlyN8M4CtAPrDOGLPdGPNswJOJY363+SiZyXFcPnWc01EkhCwry6KxvYcXd9Y6HSWqxQx3grV2N6A7UVGiqb2H1/Ye5+5lxRqbLe8zcWwKZWOT+eXacm6aV6Ab1A4ZtmhLdHluWw29HstHz6Nr5MmN6u+MZMYYPnlRKV95difrDzdywcRspyNFJRXtKHS24mqt5Werj1CQkciWyma2VDYHOZmEur+aX8B//3k/D68+oqLtEP39K+861tJFXVuXZkDKWSXEurnngmLePHCCfXVtTseJSira8q5NFY3Eug1zCjU2W87uzqUTSIpz87PVR5yOEpVUtAWAzh4P24+2MKcwg8Q4t9NxJIRlJMXx0YXjeX77MU22cYCKtgD9u9P0eixLSrOcjiJh4JMXlWCBX75V7nSUqKOiLVhr2XCkifFjEinISHQ6joSB8ZlJXDsrj6c2HaWtq9fpOFFFRVs40tBOw6lutbJlRD59cSmnuvs01DPIVLSFDUcaSYpzM6tAu9OI72YWpHPhxCx+tbacnj7tbBMsKtpRrvFUN3uOtbGoOJNYt14OMjKfvriM423dPLe9xukoUUPv0ii39nAjLmNYpq4ROQ8XT8pmRn4aP/7LIfq0j2RQqGhHsY6ePrZUNjFnfAZpiVqCVUbOGMMXrphERWMHz20/5nScqKCiHcU2lTfR67FcpOnIMgpXTs9hZkEaP1h1UK3tIFDRjlJ9Hi/rDzcyaVwKuekJTseRMGaM4b4rJlPZ2MEz29S3HWgq2lFqc2UzJ7v7WD5prNNRJAJcMW0cswrS+cGqgxpJEmAq2lGoz+PlzQMnmJCZRNnYZKfjSAQwxvClFVM42tTJExsrnY4T0VS0o9DmymZaO3u5YlqOFrIXv7l4UjYXTczm+68f1CzJAFLRjjJdvR7e2F/PhCy1ssW/jDH829VTae7o5advHHY6TsTyZY/Ih4wx5cYYa4yZGYxQEjhPb6qirauPD6mVLQEwsyCdm+YV8Mhb5Rxr0QqAgeBLS/s54GJAHVVhrrWjl/95/SCl2cmUZquVLYHxT1dNBuCbL+51OElkGrZoW2vfstYeDUYYCazvrzpIS2cv187OUytbAqZwTBKfu2wiL+yoZfWBE07HiTh+69M2xmQYY4oHH0Chv36+jM7hE6d4dF0FH1s0nrx0Lb8qgfXpi0spyU7m/ud3093ncTpORPHnjcj7gPIhxxo//nw5T9ZavvnCXhJi3fzTVVOcjiNRICHWzddvmEF5QzsPv6ltyfzJn0V7JVAy5Fjux58v5+n/dtSyal89X7hiEtkp8U7HkShx8eSxXDc7jx+sOsieY9oE2F/8VrSttS3W2orBB1Dtr58v5+fEyW7+3x93MWd8Bh+/sNjpOBJl/uMjM8lIiuMff7dd3SR+4suQv+8bY6rp759+zRizO/CxxB+stXztuV2093j47i2zidF62RJkmclxfPuvZ7Gv7iTf+/NBp+NEBF9Gj3zeWltorY2x1uZaa2cEI5iM3jNba3h5dx1f/NBkJo5LdTqORKnLp+Zw2+IiHl59WKNJ/CDG6QASGLtqWvnKsztZWprJp5aXOB1HotDgvSOn5KSSk5rA3z22hc9dNpHM5Lj3nXv7kqJgxwtb+ns5AjW393Dv41vITI7jh7fPV7eIOC4uxsUdA4X58Q2VWglwFPRujjBdvR4++8RW6tu6+cmdCzRaREJGVko8ty4az/G2Lp5+uwqP1zodKSypaEeQXo+Xv39yK+uPNPLtm2cxd3yG05FE3mdyTirXz8lnX91Jfr/lKF6rwj1S6tOOEH0eL1/87XZe21vPN26cyU3zNBlV/G9wP/X5WlqaRXevh1f2HCfO7eLGeQV+SBY9VLQjwKnuPj73xFbePHCCL189lbuWTnA6ksg5XTJlHN0eL2/sP0FHj4ebFxSSEOt2OlZYUPdImKtr7eKjP13PW4ca+K+/msXfXVLmdCQRn1w1PZdrZ+Wxt7aNO36xkYZT3U5HCgsq2mHsL/vqueb7a6hsbOeRv1nIbYs1bErCy4UTs/nY4iJ21bRy9f+s4a2DDU5HCnnqHglDnT0eHnp1P4+8Vc7U3FR+ePt8Jo5LcTqWyHmZVZDO2JR4nnq7irse2ciFE7O5Yto44mPO3F0S7WO6VbTDzFsHG/jKszupaurg7mUT+Mo109QXKGEvNz2Bz106kRd31vLWoQZ21rRy7aw8ZuSnae33IVS0w0RNSyf/9eJe/m9HLSXZyTz1qaUsK8t693l/3NUXcVJcTP9IknlFGfxx+zGe3FTFhMwkrpqRS4l2WnqXinaIa+vq5eerj/Cz1f1rEn/+ikl89tIyta4lYk3ISuZzl01kc2UTf9lXz8/XHKEkO5nlk7KZnKM1dFS0Q1RXr4dH11XwkzcP09LRy3Wz8/jyNdMoyNCuMxL53C7DkpIs5heNYWN5E2sPNfCb9ZWMS40nLsbFR+bmn7XPO9KpaIeYXo+X32+u5vuvH6SurYtJ41K4Y/EECsYk8uZ+rZAm0SXW7eKiidksK81iR3ULaw428C9/2MF3X93PbYuLuG1xETlpCU7HDCoV7RDR3t3H028f5ZE1RzjW2sX8ogyum5NHabZGhYi4XYZ5RWOYOz6DwswkfrHmCCtfO8gPVh3iymk53Ll0AheUZeFyRf5NSxVth5042c2j6yp4bEMlrZ29LC7J5Js3zeLSKWN5atNRp+OJhBRjDJdMHsslk8dS2djOkxur+N3mo7y8u46S7GRuX1zEzQsKGTNk6ddIoqIdZE9urMJrLeUN7WypbGZXTSser2V6fhq3LS6iKDOJ2tYuFWyRYUzISubL10zji1dO5uVddTy+oZJvvriXB1/dz5XTcrh+Th6XThkXcTftVbSDpLPHw5bKZp5/5xh7a9to7ewlIdbFggljuLAsm+xULaEqcj4SYt3cOK+AG+cVsK+ujac2VvHCzlpe2FlLSnwMV03P4bo5eSwrzSYxLvwLuE9F2xgzGXgUyAIagbuttRG94dtIxj2fnqFlraWjx8Pxti7q2roob2jn4PFT7EuYA+cAAAYGSURBVKxpZUd1C70eS6zbMGlcKitm5DAjP51YbVAgMiLDvTen5KYxcVwq5Q3ttHf38dKuWp7ZVkOc28W8ogwuKMtmWVkWc8dnEBcTfu8/X1vaPwV+ZK193BhzJ/AwcHngYgWGx2s51d1He3cfp04fXe///PRz26pa6O7z0N3n7T96PfR5LdaC11osAx8t/M/rB+ju83Kqq4++IQu7J8e5mZqXxicvKmVJSSaVjR1h+UIRCSdul2HiuBRuX1LEN26cybrDDaw/3Mi6w42sfP0A33utfzJPaXYyk3JSmTQuhck5KUzISmZsajxjkuJwh+hNzWGLtjFmHDAfuHLgoaeAHxpjxlprTww6LwMYuur+BIDq6uoRB/vTOzU0d/RircVr+1uxXm//517ADjzW57X0eDx0954urt7+r/u8dPV66Ozx0t7TR0ePh65ej0+/O87twu02xLldxMe4iHO7iYsxxLsMLmMwBlwDU2tdLkNxRhKxMS6S4+NJiY8hOyWO7JR4CsYkkpuWMGgabgc76mtGfC0kNFR1dwFw4pjuN4SLior+bc2K46F4eiK3TS+krTOHbUdb2FXTSnlDE+u3V/NsW9f7vs8AGUlxZCbHkRofQ0Kci8TYGBJiXSTGunG5DG5jcLsMxhhcBtymvx6YgcddxnDrovHn9df0uWqmscPsHGGMWQD8ZvAu7MaYPcCd1tqtgx57ALh/xOlERORsSqy1FYMf8OeNyJXAr4c8FgeUAh3Am8ByYOTN7shSCKxB10LX4T26Fu/Rteh3+o7pB66BL0X7KFBgjHFbaz3GGDeQP/D4u6y1LUDLGb7/gDGm+HSAof/XiDaDViyL6muh6/AeXYv36FoMb9jOFmttPbAduG3goduAbYP7s0VEJDh87R65F3jUGPP/gGbg7sBFEhGRs/GpaFtr9wFLApxFRESGEawBwy3A1zlzn3e00bXop+vwHl2L9+haDGPYIX8iIhI6NDVPRCSMqGiLiIQRvxZtY8xkY8x6Y8yBgY+TznDO14wxu40x7xhjthhjVvgzQ6jw8Vp83Bizwxiz3Riz0xjzeSeyBpIv12HQuVOMMR3GmIeCmTFYfHxNPGCMqR94TWw3xvzIiayB5uvrwhjz0YH3xq6BjznBzhpyrLV+O4BV9E9vB7gTWHWGc1YASQOfz6H/hkOiP3OEwuHjtUjjvfsKqUAlMNvp7MG+DgPPuYE3gCeBh5zO7eBr4oFI/e8/j2uxENgD5A58nQ4kOJ3d6cNvNyIHFpY6AGTZ92ZONgKT7Fkm4pj+6U8twAxrbcRMWT3Pa5EDbANWWGt3Bi9t4IzkOhhjvgp0AylAirX2S0EPHEC+XouBNXwi7r9/sBFciyeA1621v3QoakjyZ/fIeKDGWusBGPh4bODxs7kbOBxJBXuAz9fCGHODMWY3/a3sByOlYA/w6ToYY2bT/xfY94KeMHhG8v742EC32avGmGXBDBkkvl6L6UCpMWa1MWarMebfzaB57tHKsRuRxphLgG/w3vT4qGStfd72r6A4GbjLGDPF6UzBZIyJBX4O3Hv6TRzlfkr/ym6zgQeBPxpjshzO5JQYYDb9y0JfAlwN3OVoohDgz6L97sJSAGdbWGrguWXA48CN1tr9fswQKny+FqdZa6uATcB1QUkYHL5chzygDHjRGFMB3Ad8yhjzsyBnDTSfXhPW2jprbe/A538eeH5mkLMGmq/vj0rgD9babmvtSeCPwOKgJg1Bfiva1seFpYwxi4DfAjfbQetxR5IRXIupgz7PBi4DIqZ7xJfrYK2tstZmW2uLrbXF9C/x+3Nr7aeDHjiARvCaKBj0+VygGIioho2v14L+m9JXmX6xwBXAO8FLGqL8eVcTmApspP8mw0ZgysDjLwILBz5/GzhB/z/a6WOW03dk/X34eC2+B+weuAbvAP/gdG4nrsOQ8x8gQkdP+PiaeBTYNfB6eBu4xuncDl4LF/DfwN6B98l/Ay6nszt9aBq7iEgY0YxIEZEwoqItIhJGVLRFRMKIiraISBhR0RYRCSMq2iIiYURFW0QkjKhoi4iEkf8Pua19dAfCeGUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p_thresh = [(sample>np.log(4)).mean() for sample in pp_samples['dist']]\n",
    "ax = sns.distplot(p_thresh)\n",
    "ax.vlines((hennepin_radon>np.log(4)).mean(), *ax.get_ylim(), color='r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shows a histogram of the proportion greater than log(4) calculated for each simulated dataset, compared to the empirical proportion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prior sensitivity\n",
    "\n",
    "Its also important to check the sensitivity of your choice of priors to the resulting inference.\n",
    "\n",
    "Here is the same model, but with drastically different (though still uninformative) priors specified:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 2 jobs)\n",
      "NUTS: [σ, μ]\n",
      "Sampling 2 chains: 100%|██████████| 3000/3000 [00:02<00:00, 1396.83draws/s]\n"
     ]
    }
   ],
   "source": [
    "from pymc3 import Flat, HalfCauchy\n",
    "\n",
    "with Model() as prior_sensitivity:\n",
    "    \n",
    "    μ = Flat('μ')\n",
    "    σ = HalfCauchy('σ', 5)\n",
    "    \n",
    "    dist = Normal('dist', mu=μ, sd=σ, observed=hennepin_radon)\n",
    "    \n",
    "    sensitivity_samples = sample(1000, cores=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3wUdf7H8dfsphdSSEIqKRB671WqKKByAmIDRM4u6FnO49CfoN7ZUE8QO4qi2D0QFQSkI72X0EJCekhCSM+m7M7vjwntJCSQMls+z8djHslmZ2c/+4Xd936mKqqqIoQQQtgbg94FCCGEEA1BAk4IIYRdkoATQghhlyTghBBC2CUJOCGEEHZJAk4IIYRdkoATQghhlyTghBBC2CUJOCGEEHZJAk6IBqYoiqooSsuLbn+mKMq/9KxJCEcgASeEEMIuScAJIYSwSxJwQggh7JIEnBCNw+2i3311q0IIByIBJ0TjuFdRFKOiKF2AYYC3oijOehclhD2TgBOicXgAGcDHwPPAPcBQXSsSws4pcsFTIRqWoigqEKuqarzetQjhSKSDE0IIYZck4IQQQtglWUUphBDCLkkHJ4QQwi5JwAkhhLBLTjXcL+svhRBCWDOlujukgxNCCGGXJOCEEELYJQk4IYQQdkkCTgghhF2SgBNCCGGXJOCEEELYJQk4IYQQdkkCTgghhF2SgBNCCGGXajqTiRCiBpVmC8dPF3E4PZ/0PBPZRSYqzSpORoXm/h7ENvOmZ5Q/Xq7ydhOiMck7TohroKoqO0+d5btdKfx2KJOissrz9/l5OONsNFBWaSG/tAIAo0Ghe6Qf47uFM7pTCJ4SdkI0uJoulyPnohTiIqqqsvxgJu9viOdQWgFerk6M6hhM3xZN6RTuS7ifO65OxvPz55WUE5dewB8nc1hxKJOE7GKauDnx4KAW3Ns/Cg8XCToh6qjac1FKwAlRSydOF/Ls0kPsSMwlJsCT+6+LYUyX0FqHlKqq7E46y/vrT7LmaBZB3q48f3M7RncMQVGqfY8KIa5MAk6Ia2WxqHz6RyKv/XYUT1cnZtzYhtt6RGA0XHso7U7KZdaywxxKK2B422a8OaEzPu7O9Vi1EA5DAk6Ia1FoquCxr/ey7lg2I9o145WxHWnq5VovyzZbVBZWBWeYrzsfTe5Bq2be9bJsIRyIBJwQVyu7sIwpC3dwLLOQ529ux6Q+kQ2yKnHXqVweXrwHU7mZ9yd2Z0BsQL0/hxB2TAJOiKuRklvCpE+2k1lg4v2J3RnSOqhBny89r5Spn+0kPquI18d3Ymy38AZ9PiHsiAScELUVl17APQt3UF5p4dMpPeke6dcoz1tgquChL3azNeEMr43txISeEY3yvELYOAk4IWrjSEYBd3y0DXdnI4v+2qvRt4mZKszcv2gXm07k8Pr4TkzoISEnRA0k4ISoyamcYsZ/sBUng8L3D/Ulwt9DlzrOhdyWk2f45J4eDG7g1aNC2LhqA07ORSkEkFtczqRPt2NRVb68r5du4Qbg5mzk/YndadXMm0cX7yEuvUC3WoSwZRJwwuGVV1p46MvdZBWU8emUnrQM0n9XfS9XJxZO6Ym3mzMPfLGLvJJyvUsSwuZIwAmHN2uZdnaS18d3okuEr97lnBfs48YHk7qTVVDG49/sw2KRLQZCXA0JOOHQluxN5esdKTw8uAVjuoTpXc6fdInw5fmb27HheDbz18XrXY4QNkUCTjishOwinl1yiJ5Rfjx1fSu9y6nW3b2bM6ZLKHPXnGBP8lm9yxHCZkjACYdkqjAz7au9uDoZmHdnV5yM1vtWUBSFF8d0ILiJG098u++SS/MIIapnve9qIRrQy8uPEJdRwJsTOhPi4653OTXycXfmP7d3ISW3hJeXH9G7HCFsggSccDi/Hcpg0dYk7h8YzdA2zfQup9Z6Rfvz1wHRfLU9ma0nz+hdjhBWTw70Fg4lu7CMEf/ZQHN/D75/qB8uTrb1Ha+03MzIuRtRgd8evw53F2ONjxHCzsmB3kKAdkhAcbmZNyd0trlwA3B3MfLquE4knSnhnbUn9C5HCKtme+9wIa7RrwcyWH4wkyeGt7KKg7mvVZ+YpoztGsaCzYmk5JboXY4QVksCTjiEQlMFs38+TMcwH+4fGK13OXX29xtbY1QUXl1xVO9ShLBaEnDCIfxn9Qlyisr4960drPqQgNoK8XHnoUEt+PVgBjsSc/UuRwirZPvvdCFqcCSjgM+3nuKuXs3pFG49p+KqqweuiyHEx40Xfzksp/ES4jIk4IRdU1WV5386hI+7M3+/obXe5dQrdxcj/7ixDYfSCvhxT6re5QhhdSTghF377540dp46y4wb2+Dr4aJ3OfXuls6hdInw5fWVxyiWM5wIcQkJOGG38ksreGXFEbo292V893C9y2kQBoPC8ze3I7uwjE83J+pdjhBWRQJO2K3/rD5ObnE5L43pgMFQ7bGgNq9bcz+Gtw1iweZECkwVepcjhNWQgBN2KTGnmC+3JXF7zwg6hPnoXc6frFu3jiFDhuDj40NUVNQV542Li6NHjx74+fnh5+fH8OHDiYuLu2Sex4e1Ir+0gs//ONVwRQthYyTghF16/bejuDgZeMJKL4Pj6enJ1KlTmTNnTo3zhoaG8sMPP5Cbm0tOTg633HILd9xxxyXzdAz3kS5OiP8hASfsRlRUFHPmzCG2TXs+nNqfJjs/QS3JZ+TIkXh7ezN8+HDOntWup7Zt2zb69euHr68vnTt3Zv369eeXs3DhQtq2bYu3tzcxMTF8+OGH5+9bv3494eHhvPnmmwQFBRESEsLChQuvutZevXoxadIkYmJiapzX19eXqKgoFEVBVVWMRiPx8X+++Kl0cUJcSgJO2JUff/yRFve8SucnPuPk7g2MHDmSl19+mZycHCwWC/PmzSMtLY3Ro0fz3HPPkZubyxtvvMG4cePIzs4GICgoiF9++YWCggIWLlzIE088wZ49e84/R2ZmJvn5+aSlpfHJJ5/w6KOPng/OV199FV9f32qnuvD19cXNzY3p06czc+bMP90vXZwQl5KAE3ZlwJiJHM038M/b+nPdwIH07t2brl274urqyq233srevXv58ssvGTVqFKNGjcJgMHD99dfTo0cPli9fDsDo0aNp0aIFiqIwaNAgRowYwaZNm84/h7OzM88//zzOzs6MGjUKLy8vjh07BsCMGTPIy8urdqqLvLw88vPzmT9/Pl27dr3sPNLFCXGBBJywGyqwMsFEm2BvxnULx93dnWbNLlzvzd3dnaKiIpKSkvj+++8v6aw2b95MRkYGACtWrKBPnz74+/vj6+vL8uXLycnJOb+cpk2b4uTkdP62h4cHRUVFjfIaPT09eeihh5g8eTJZWVl/uv9cF/fxpgTp4oTDk4ATdqO4rJKswjKeHd0W4xUOC4iIiGDSpEmXdFbFxcXMmDGDsrIyxo0bx9NPP83p06fJy8tj1KhR1HDdxPNefvllvLy8qp3qg8VioaSkhLS0tMve/9iwWApMlXyzI7lenk8IWyUBJ+xCcVklhaZK2od4MzA28IrzTpw4kZ9//pmVK1diNpsxmUysX7+e1NRUysvLKSsrIzAwECcnJ1asWMGqVatqXcfMmTMpKiqqdjrHYrFgMpmoqKhAVVVMJhPl5eWXXebq1avZu3cvZrOZgoICnnzySfz8/Gjbtu1l5+8U7kufGH8W/nGKCrOl1rULYW8k4IRd+HzrKcwWlfE9ImqcNyIigp9++omXX36ZwMBAIiIimDNnDhaLBW9vb+bNm8eECRPw8/Pjq6++4pZbbqn3ejdu3Ii7uzujRo0iOTkZd3d3RowYcf7+9u3bs3jxYkDb9nbnnXfi4+NDixYtiI+P57fffsPNza3a5T9wXQwZ+SZ+PZBR77ULYSuUGla9yCnKhdUrMFUw8LV1dI/049MpPfUuxypYLCoj3t6Iq5OBX6YPQFHs90wuwuFV+59bOjhh8z7ZlEh+aQVPWulB3XowGBTuGxDN4fQCtp48o3c5QuhCAk7YtLPF5XyyOZGRHYKt8pRcevpL1zACvFz4eFOC3qUIoQsJOGHTPtyYQHF5pdWekktPbs5GJveNYt2xbE6cLtS7HCEanQScsFlZhSY+25LImM6htGrmrXc5Vmlin0jcnA0s2CSX0hGORwJO2Kz315+kwqzy+HDp3qrj7+nC+O7hLNmbRnZhmd7lCNGoJOCETUrPK2XxtmTGdwsnOsBT73Ks2r39oyk3W+TAb+FwJOCETZq/Tjub/mPDY3WuxPq1CPRiYGwAX25PkgO/hUORgBM2J/lMCd/tTOHOXhGE+brrXY5NmNIvitMFZaw6fFrvUoRoNBJwwua8tz4eg0HhkSEt9S7FZgxuHUSEvzufbz2ldylCNBoJOGFT0vJK+XFPKnf0jKBZk+pPVSUuZTQoTO4TxY7EXI5kFOhdjhCNQgJO2JQPN5wE4MFBLXSuxPbc1iMcN2cDi7ae0rsUIRqFBJywGVkFJr7ZmcK4buGy7e0a+Hq4cGvXMJbsTSOv5PJXLhDCnkjACZvx0cYEzBaVRwbLtrdrNblvFKYKC9/vSr32hez7WpuEsHIScMImnCkqY/H2ZMZ0CaV5Uw+9y7FZbUOa0Cvan0XbtMsLXZP8FG0SwspJwAmbsGBzIqZKs3Rv9WBy30hSckvZeDxb71KEaFAScMLq5ZWUs2jLKUZ3DKFlkJfe5di8Ee2CCfR25YttSXqXIkSDkoATVm/hH6coLjczbah0b/XBxcnAnT0jWHcsi5TcEr3LEaLBSMAJq1ZoqmDhH4mMaNeMNsFN9C7HbtzZuzkGReErOT+lsGMScMKqLdqaRIGpkulD5ZyT9SnEx53hbYP4dmcKZZVmvcsRokFIwAmrVVJeyYJNCQxpHUjHcLlad32b2CeS3OJyVhzM1LsUIRqEBJywWou3JXO2pIJp0r01iP4tAogO8JSdTYTdkoATVslUYebDjQn0b9mU7pF+epdjlwwGhbt7N2d30lni0uX8lML+SMAJq/TtzhRyispk21sDu617BG7OBr7cLl2csD8ScMLqlFWa+WDDSXpF+dMnpqne5dg1Hw9nbukcytK9aRSYKvQuR4h6JQEnrM6Pu9PIyDfJcW+NZGKfSErKzSzZk6Z3KULUKwk4YVUqzBbeWx9P5whfBsYG6F2OQ+gU7kvncB++2JaEql7j+SmFsEIScMKq/LQvndSzpTw2tCWKouhdjsOY2CeS+KwitiXk6l2KEPVGAk5YDbNF5b118bQLacLQNkF6l+NQbu4cio+7s+xsIuyKBJywGr8ezCAhp5jp0r01OjdnI7d1D2floUyyCkx6lyNEvZCAE1bBYlGZv/YEsUFe3NA+WO9yHNLdfSKptKh8s1Ou9SbsgwScsAqr4jI5frqIaUNbYjBI96aH6ABPBsYG8NX2ZCrNFr3LEaLOJOCE7lRV5Z218UQHeHJTp1C9y3Fok/pEkllg4vcjWXqXIkSdScAJ3a07lsXh9AIeGdwCo3RvuhraJohQHze+lPNTCjsgASd0paoq89bEE+7nzl+6huldjsNzMhq4s1dzNsfnkJBdpHc5QtSJBJzQ1eb4HPal5PHw4BY4G+W/ozW4vVcETgaFxdvlYqjCtsknitDVO2vjCW7ixvju4XqXIqoEebtxY4dgvt+VQmm5XAxV2C4JOKGb7Qln2JGYy4ODYnB1MupdjrjIpD6RFJgq+Xl/ut6lCHHNnPQuQOgo4wD8PgvS90KFCXybQ6/7teliJbnwXl8oygRXH/jnFVZdHfwBdi6AnONQXgJ+UdD3Ueg2Sbs/NxGWPASZB/BVYujo+Sh39mqu3bfuZTi+Eu5fB4Z6+O71+2w48jOcidduj3kPut595ccsGgOZh8CUD+6+0LwPjPiX9jpAG4uVz8LJtVCaC55B0PYmuP5FcHLV5qkohXX/hkNLoOg0eDSFrhNh2P9p9y+bDsnbIT8VnFwgrAeMeAmC2lY93gQbX4eD30PhafCPgSH/hHZjGn5M1r0CG16lN3DKDfi1avp7AnhWc2WHmsYkfS/8NB1yEyBqANz6AXj4a49d8hCUFcIdi+v+2oT4H9LBObJv7tY+lLxDIfZ6LZSWPw2JGy+d75e/QUlO7ZYZvwbOnoIWw7RwyD4Cy6bBsRXa/av/DzIPkBs+lOZlx5nr/yNuzkYt+P6YB6PfvPZwK8mF1F0XbqfuAp8IcPev/TIK0qHlcO1D3+CkhcHSRy7cv/JZ2P8VoEK7v2gfzts/gE1vaverKnw7Eba8A0Yn6HKnNg65CReWsWcRuHpDx3Haz/jV8MVYLdgAVs7Ulmdwhs53aF8svrsHUnbUov4MyDxY/f21HZO2txDX/G4+rbyR7PZTwdmt+nlrGpOfH4eCVO3/WPxq2PyW9vfkbRD3E9z4Ss2vS4hrIB2cozJXaB86AOMWQLN28OEgyNgHeRd1aPu+0j7kr3sGNrxa83J73Q83z9U6E4CFoyFpM5xcB61HQvZxiL6Op8v/xr1KGv0tVc/12wzoOB7Ce1zd67CYtVDdt1gL0Q7jLixjyi/az/cHaJ1FbUzbeeH3uGXw3SQ4e9Eu82cTtZ8Dn4beD2gf7lvnXxizxA0Q/zsEtIIHN10+GKaugua9q5aXBHM7QWE6ZB+F0C4Qt1S775Z5WscT2AZW/lMLjLu+/fPyKsvg6K/av9XJtXDd3yG44+VfX23HpNcDRIT0YfzLa4gjhDdcPKuft6YxyT6udfCj5sC7vSHrqPbv9uvTMOBJbc2BEA1AAs5RGZ2h98Ow7V347/3aarCM/dCsI7S5SZsnLxlW/AP6TtM+aDfUYrlh3S69bS7XfjapOoA7sBXmE2sYU1ZAH+cjGJrdAMd+077NT99T+/qzj2uhtv8brcMJageDZ0Cn22u/jOqsf00LnOOrQDFC/8cu3Nf7Ia0L2vQGpO7UVqm6+19YrZtQNUgunvBBfyjMhJAuMPLVC6FzLtxA+6IBoBjAu+oUZU5VoZixH8K6w+lD2u3MQ5fWmbZHG4ODP4ApT5v3xleg4211H4Nv7sbbXM4aj1DeODCSs6Pa4ufpcvl5axqTwFZw6EcoztHWEsRer63Grii+dGyFqGcScI6szWg4+rP2AXr6kLZKrM1obbWZxaJtH/GNhKH/Bynbr375W+ZD6g4tPHtM1f52/UskJiYyonI3SmhXGPIsfDVB2z51/DfY+q42X99HtO1Wl7N4ApxYCd4hWtfX6XYI6XRtY3A5e7+E/KruI6A1hHS+cF9Yd4jorXWlB7/T/tZuDPhFa7+XnNF+pu+F1qO1Gk9t0mqevksLvnPKimDpw1Wvd9qFgBv4JPz6lLaqcuXMC/MXnb7w+7tVq3/9oqH3g9oYNG1R99ducILI/hAQC3nJhJxcy5vG+Sxf1YZRt066/GNqGpOb52rb4I6vhJbXQ9fJ8MlwGPeJ9u+9/2ttXAbPhFYj6v4ahKgiAeeoSnJh8XioKIF7f4OgNtp2oA2vgmcAtLoBkv6AZh20bUrnVmdVFGsf1mPeBa/A6pe//lVY/4q2c8bkZeDWBIBj5QHckDeTx4bF8uT1rbSdGtx8tA/V9/tp2+BUi7YjRngv7dv//8qvWrUa2Fqrzz+6fsfmiYNQXqx98P76lBbATx7RPoS/vwfSdsOwWdDnYVj1nNaNoMCEz7WxA2214p1fgbkSXo/WOsL0fRDVX7u/+Iw2/ul7oNs92g4Z5/S8D0K6aqsbUbXgWzb9wrLPj4ECzdprY+BTT4dZXPc0DPr7hds/TIVDP1JxeBnmMRMvf6aZmsYktCs8vPnC/EsehsgBWte65gW463tI2QbfT4EnD4O7X/28FuHwZCcTR3X2lBZuBmdttaK7nxYYoK1GOndl59OHtG4ptWrblKVSu11Rot3OPq5NlWVV91vglye1cAvupG1v8o04/7Tz18Xj6WJkav8orYYt82DUm5B1RAu2sO7apFog6/Dla39wA9z+JTi5azuwzImF7+/VtsGdW+VXG8VntNoLqnaFLyu68LpdPC+sqjXlQ3G29nvWUe1neA9wdtc+vM+NGWiBU51z3VteMnw6Qgu3AU9o29ouvjxQZTmEd9eCZtAzkLRV+3vM4AvzPH1c2xvRlA/fTdbG4KdHtVWkljqcKPninWEuUmCy8PuR09r4Fudo07nnqWlMLpa8XdvGeOMrF3aGieipdckVxdU+vxDXQjo4RxXYWgu10rPw+S1aF3TwB+2+5n3ALxJm51+YP3ETfH7Tnw8TeLen9vPBTdpqwnX/gl2faN/OQzrB5v9o9/vHcDLmLn45kM6D17XA18MFls6ADmO1D7hzO2P8+jRQFTJNYy9fu9EZ2t6sTUVZcOBbbQeLr++ArpNgzHxtvk1vQc4JyK+6/MueRXBqM3SbDJF9YcdHWsfaerTWbcX9pO2eH95L+6BOWHehDt/IqrHprXVWy6ZrgXP01wtjBtC2atVc9lH4+i4oL4SyAm37W7MO2jyfjIDCDG1vxgoTrJih/b3jbVqw7fkcDnyn7fiTdURbPezqo+08co6Lh7aHZec7tB1V9n+tTXu/hCHPXdqFXaymMfniL+DVTNummZ8KJ9egKga2eQwiZ3MiN4SWwY4PtccOmVl1KEUNY3KOxQzLn9JC3S9S2xEHtC6xIB2MLuAbdfm6hbgGEnCOysUT7v4B1r6k7cyQsb9qW9m92p6I16ogQ/upWrQP23MiB/BeUm9cnQzcNzBa2x6TvOXCjiXBHWH4bG27HWi/B3eo+fm8gqDfdG1K2wMFaRfui1+jbRc6J2WbNkUN0D7M/1fTluAZWNWhmrRld50Ig2Zc6LD+8j6sngUJ67VQ9QyCnvdr9YJ2aMDEH7WdcxLWgbOHNp4j/qXdB1q4gRYy29+/8PzBHbWA84/Wvnjs+1r70G91o7b86rax+UVqO9gM+ocWVpbK6serpjHpdg8cXqLtFKIYIaI3ysCn6JzZkpeXH+V4ljd/Wmlc05ics+tTrUvu/7h2u/VIbUen/V+BsyfcPK/6Y+2EuAaKem6VzOVd8U4haiv5TAlD3lzPPX2jeP7mdnqXI65SfkkFfV9dw6iOIbwRtFL746Bn9C1KCE21lyCRbXCiUby/IR6jQeHBQTF6lyKugY+HM+O6hbNsXzrFZVfoEIWwIhJwosGl55Xyw+5Ubu8RQbMmVzgjhrBqU/pHUW62cCA1v+aZhbACEnCiwX244SSqinRvNq5FoBdDWgeyLzWPyrrsqSlEI5GAEw0qq8DE1ztTGNctnHA/D73LEXU0dUA0JWWVHMss1LsUIWokASca1MebEqg0W3hkSD2cZUPobkDLAAK8XNiddJYadlATQncScKLBnCkq48ttyYzpEkZk0yucrFfYDEVR6BHlT3ZhGeuPZetdjhBXJAEnGsyCzYmYKs08OqSl3qWIetQmuAnebk68v/6k3qUIcUUScKJB5JWUs2jLKUZ3DKFlkJfe5Yh6ZDQodI/0Z8epXHYn1fIyRELoQAJONIhP/zhFcbmZaUOle7NHHcN98PVw5v31cu5IYb0k4ES9KzBVsPCPRG5o34w2wU30Lkc0ABejgXv6RvH7kdOcOC17VArrJAEn6t2iLacoNFUyfWg1J0sWduGeflG4ORv4YIN0ccI6ScCJelVcVsknmxMZ2iaIDmE+epcjGpC/pwt39GzOT/vSSM8r1bscIf5EAk7Uqy+3JXG2pILpsu3NIdw3MBoV+GRzot6lCPEnEnCi3pSWm/l4UwIDYwPo2lyuyuwIwv08GNM5lK93JJNXUq53OUJcQgJO1JuvdySTU1Qu294czIODWlBSbmbR1iS9SxHiEhJwol6YKsx8uPEkvaP96RXtr3c5ohG1DvZmWJsgPttyitJys97lCHGeBJyoF9/vTuV0QRmPDZPuzRE9NLgFucXlfLMzWe9ShDhPAk7UWXmlhQ/Wn6Rbc1/6tWiqdzlCBz2jtM79gw0nMVVIFyesgwScqLP/7kklLa+U6cNiUZRqrx4v7Nzjw2I5XVDG97tT9S5FCEACTtRRpdnCe+tP0inch8GtAvUuR+ioX4umdI/04/118ZRXygVRhf4k4ESd/LQvneTcEqYPle7N0SmKwmPDYknPN/HjHunihP4k4MQ1M1tU3l0XT9uQJgxvG6R3OcIKXBcbQOdwH95bH0+FWbo4oS8JOHHNfj2YQUJOMdOHtpTuTQAXuriU3FKW7k3Tuxzh4CTgxDWxWFTmrz1BbJAXN7YP1rscYUWGtgmifWgT3l0XT6V0cUJHEnDimqyKy+T46SKmDW2JwSDdm7jgXBd36kwJPx9I17sc4cAk4MRVU1WVd9bGEx3gyU2dQvUuR1ih69s2o02wN++sjcdsUfUuRzgoCThx1dYezeJwegGPDG6BUbo3cRkGg9bFJWQX8+vBDL3LEQ5KAk5cFVVVmbc2nnA/d/7SNUzvcoQVu7F9MLFBXsxfewKLdHFCBxJw4qpsOpHD/pQ8HhncEmej/PcR1TMYFKYNbcnx00WsOJSpdznCAcknlKg1VVWZt+YEIT5ujOsu3Zuo2U2dQokJ9GTeGuniROOTgBO1tjXhDLuSzvLQoBa4Ohn1LkfYAKNB4fFhsRw7XcjKw9LFicYlASdq7Z018QR6u3J7zwi9SxE25FwXN1e6ONHIJOBErew6lcvWhDM8eF0Mbs7SvYnaMxoUHhsay9FM6eJE45KAE7Uyb208/p4u3NW7ud6lCBt0c+dQYgKkixONSwJO1GhfSh4bj2dz/8AYPFyc9C5H2CBj1XFxRzMLWRUnXZxoHBJwokbz157A18OZSX0j9S5F2LBzXdzbv0sXJxqHBJy4osPp+fx+JIup/aPxcpXuTVw7o0Fh+rCW0sWJRiMBJ65o/tp4vF2duKdflN6lCDtwc6dz2+LipYsTDU4CTlTr+OlCVhzKZEr/KHzcnfUuR9gBJ6OB6cNaciSjgFVxp/UuR9g5CThRrflr4/F0MTK1f7TepQg7cnOnUKJlj0rRCCTgxGWdzC7ilwPpTOwbiZ+ni97lCDviZDQwfah0caLhScCJy3r79xO4ORu5f2CM3qUIO3RLZ62Lm7fmBKoqXZxoGBJw4k+OZBTw8/507u0fRYCXq97lCDvkZDQwbUhL4rwQ7WEAABtuSURBVKSLEw1IAk78yVurj+Pt5sQDA1voXYqwY2O6hBLV1IO5v0sXJxqGBJy4xP6UPFbHneb+gTH4eMiek6LhaNviYonLKGC1dHGiAUjAiUu8ufo4fh7O3Ns/Su9ShAM418W9LV2caAAScOK8HYm5bDyezUODWuDtJt2baHhORgPTpIsTDUQCTgDa1brfWHWMQG9XJveN0rsc4UD+0iWUyKYezJU9KkU9k4ATAPwRf4YdiblMG9ISdxe53ptoPOe2xR1OL+D3I1l6lyPsiAScQFVV5qw6RqiPG3f0kqt1i8Z3rot7+/fj0sWJeiMBJ1hzJIv9KXk8NiwWVyfp3kTjO3dcnHRxoj5JwDk4i0XlzdXHiWzqwbju4XqXIxzYrV3DqrbFSRcn6ocEnIP7+UA6RzIK+NvwWJyN8t9B6OdcF3corYA10sWJeiCfaA6srNLMnJXHaBvShDGdw/QuRwhu7RpGc38P3pYuTtQDCTgHtnhbMqlnS5kxsg0Gg6J3OUJUHRcnXZyoHxJwDqrAVME7a0/Qv2VTrosN0LscIc4718XJcXGiriTgHNRHGxI4W1LBjBvboijSvQnr4Vy1Le5gWj5rj0oXJ66dBJwDyiowsWBzAjd3DqVjuI/e5QjxJ7d2CyPC313OUSnqRALOAb295gRmi8rfR7TWuxQhLsvZaGD6kFgOpuWz7ph0ceLaSMA5mGOZhXyzI5m7e0fSvKmH3uUIUS3p4kRdScA5EFVV+devcXi7OfP4sFi9yxHiis5tizuQKl2cuDYScA5k7dEsNp3I4W/DY/HzdNG7HCFqNLZbOOF+0sWJayMB5yDKKy38+9cjtAj0ZGKfSL3LqZMjR44wdOhQfHx8aNmyJUuWLLnsfC+88AKKovD777+f/9tXX31FSEgI0dHRrF+//vzfT548Sb9+/TCbzdU+72effcaAAQP+9PeoqKjzz/HZZ59hNBrx8vKiSZMmdOnShV9++QWA9evXYzAY8PLywsvLi/DwcCZMmMDOnTuvZRgcgrPRwPShWhe3/li23uUIGyMB5yAWbT1FQk4xz93UzqZPyVVZWcmYMWO46aabyM3N5aOPPmLixIkcP378kvlOnjzJDz/8QEhIyCWPnTFjBnv27OGdd95h2rRp5+977LHHeOuttzAa636y6b59+1JUVEReXh5//etfmTBhArm5uQCEhoZSVFREYWEh27Zto02bNgwcOJA1a9bU+Xnt1YUuTs5uIq6O7X7SiVrLLS5n7poTDGoVyJDWQXqXUydHjx4lPT2dJ554AqPRyNChQ+nfvz9ffPHFJfNNmzaN1157DReXC6tiz5w5Q1hYGCEhIQwfPpyEhAQAfvjhB8LCwujTp0+91mowGJg6dSqlpaXnn+scRVEIDw/nxRdf5L777uMf//hHvT63PTm3LW6/dHHiKknAOYC3Vh+jpNzMc6Pb6l1KnV3uG7yqqhw6dOj87e+//x4XFxdGjRp1yXyBgYGcOXOG1NRUVq9eTfv27SkqKuJf//oXr7zySr3XWllZyYIFC/Dy8iI2tvqdesaOHcuePXsoLi6u9xrsxfkuTs5uIq6CBJydO5Cax+LtyUzqE0lsM2+9y6mzNm3aEBQUxJw5c6ioqGDVqlVs2LCBkpISAIqKipg5cyZvv/32nx5rMBh4//33GT9+PG+88QYff/wxzz//PNOnT+fgwYMMGTKEG2644ZKw/F/btm3D19f3kik5Ofmy8wQHB/P111+zZMkSfHyqP6A+NDQUVVXJy8u7xlGxfy5OVV1cSh7rj0sXJ2rHSe8CRMMxW1SeXXKIAC9XnhzRSu9y6oWzszNLly5l+vTpvPbaa/To0YMJEybg6uoKwKxZs5g0aRLR0dGXffywYcMYNmwYAAcOHGDXrl3MmTOHqKgoNm/eTEpKCvfddx/btm277OP79OnD5s2bL/lbVFRUjfNcSVpaGoqi4OvrW+vHOKKx3cJ5Z208b/9+gsGtAuUUc6JG0sHZsS+3JXEwLZ/nb2pHEzdnvcupN506dWLDhg2cOXOGlStXkpCQQK9evQBYs2YN8+bNIzg4mODgYFJSUpgwYQKvvfbaJctQVZVp06Yxb948cnJyMJvNREZG0rNnTw4cONCor2fJkiV069YNT0/PRn1eW+PipF1pQLo4UVvSwdmprAITb6w8xsDYAG7qFFLzA2zIgQMHaNWqFRaLhffee4+MjAymTJkCaAFXUVFxft6ePXvy1ltvMXLkyEuWsWDBArp27UqXLl2orKyktLSUuLg4kpOTiYmJafDXoKoq6enpLFiwgAULFrBs2bIGf057MK5bOPPXxjNXujhRCxJwduqlX49QZrbw4pgOdvch8MUXX7BgwQIqKioYOHAgq1evPr+KsmnTppfMazQa8fPzw8vL6/zfcnJymDt3Llu2bAHAycmJ+fPnM3ToUNzc3Fi4cGGD1Z6eno6XlxeqquLj40O/fv1Yv359ve/Baa9cnAw8OqQlM5ccZMPxbAbb+F7BomEpNeyRJLsr2aCNx7OZ/OkOnhjeiseHyym5RD3b8Lr2c9Azujx9eaWFIW+sJ9DblSWP9LO7L3DiqlX7H0C2wdmZ4rJKnl16kOgATx4a3PCr2oRobOe6uH0peWw8kaN3OcKKScDZmVdWHCH1bCmvj++Eq1Pdz8ohhDUa3z2cMF85u4m4Mgk4O7L5RA5fbkvmr/2j6Rnlr3c5QjQYFycDjwxpwd5k6eJE9STg7EShqYJ//HiAmEBPnr5BLmQq7N9t3SMI83VnzsqjWCzSxYk/k4CzE//+9QgZ+aW8eVtn3Jxl1aSwfy5OBp6+oRWH0gr4aX+a3uUIKyQBZwfWHcvim50pPDioBV2b++ldjhCNZkznMDqENWHOb8cwVVR/qSPhmCTgbFxWgYmnv9tP62be/E0OCRAOxmBQmDmqLen5Jj79I1HvcoSVkYCzYWaLyuPf7KOk3Mz8u7rKXpPCIfVrEcDwtkG8t+4kZ4rK9C5HWBEJOBs2b80Jtiac4cUx7e3iSgFCXKsZI9tQWmFm3poTepcirIgEnI3aEp/DvLUnGNstjNt6ROhdjhC6ahnkzZ29Ili8PZmT2UV6lyOshAScDcouLOPxb/cRE+DJS2M66F2OEFbhb8Nb4eZs5MWf4+TgbwFIwNmc8koLjy7eQ0FpBe/e3Q1PVzlfthCAdt3D61ux4Xg2Kw9n6l2OsAIScDZEVVWe/+kQO07l8vr4TrQJbqJ3SUJYlcl9I2kb0oQXfo6juKxS73KEziTgbMhnW07xzc4UHh3SgjFdwvQuRwir42Q08K+/tCcj38S8tbLDiaOTgLMRm05k89IvcVzfrhlPXS+n4hKiOt0j/ZnQI5xPNiVy/HSh3uUIHUnA2YCT2UU8ungPrZp58/btXTAY5PpXQlzJjJFt8XJz4rmlh2SHEwcmAWflMvJLmfzJDlycDHw8uYfsVCJELfh7uvCPG9uwIzGXH3an6l2O0IkEnBU7W1zOpE92UFBawWf39iLC30PvkoSwGbf3iKBnlB8v/hJHZr5J73KEDiTgrFRxWSVTPttJSm4JC+7pQYcwH71LEsKmGAwKc8Z3psJsYeaSg7Kq0gFJwFmhskozD325m0Np+cy/qxu9Y5rqXZIQNikqwJNnbmjD2qPaFTeEY5GAszKmCjMPfbGbTSdyeG1cJ65v10zvkoSwaVP6RTEwNoDZyw5zNLNA73JEI5KAsyKmCq1zW3csm5dv7cj47uF6lySEzTMYFN6a0IUm7s48uniPHADuQCTgrISpwswDX+xm/bFsXh3bkbt6N9e7JCHsRqC3K3Nv70JCTjH/J4cOOAwJOCtQWm7m/kW72HQim9fHdeKOXhJuQtS3fi0DeHxYLP/dm8b3cuiAQ5CA09nZ4nLuXrCNzfHaNrcJPeXSN0I0lOlDY+nXoinP/3SIuHTZHmfvJOB0lJ5Xym0fbuVQWgHv3dWNCXJdNyEalNGg8PYdXfB1d+G+z3eSVSjHx9kzCTidHD9dyLj3t3A638Siv/ZiZMcQvUsSwiEEebux4J4enC2p4MEvdlNabta7JNFAJOB0sDspl9s+2EqlReXbB/vSR45zE6JRdQjz4T+3d2FfSh4PfrmbskoJOXskAdfIfo87zV0fb8ff04X/PtyPdqFyTTch9HBjh2BeHduRjcezmf7VXirMFr1LEvVMAq4RfbEtiQe+2EXrYG9+eKivnFtSCJ3d3rM5s29ux6q40zz13X7MFjl8wJ7IqekbgcWi8upvR/loYwJD2wTxzp1d5aoAQliJKf2jKakw8/pvx3B3NvLK2I5ySSo7IZ+yDcxUYebJ7/ax/GAmk/pEMuvmdjgZpXEWwpo8MrglpnIz89bGU1ReyVsTOuPqZNS7LFFHEnAN6ExRGfct2sW+lDyeG92Wvw6IRlHkm6EQ1ujJEa3xdnPm38uPcKaojI8m96CJm7PeZYk6kFaigZzMLuLW97YQl17A+3d3476BMRJuQli5+6+L4e3bu7Dr1Flu/3AbWQVynJwtk4BrAFtO5jD2vS0Ul1XyzQN9uLGDHOMmhK34S9cwPp3Sk6Qzxdz63hYOpeXrXZK4RhJw9eybHclM/mQHQd6uLH20P12b++ldkhDiKl3XKpDvHuyLRVUZ/8EWftqXpndJ4hpIwNUTs0Xl5eVHmPHfg/RrGcCPj/STwwCEsGEdwnz4efoAOoX78vg3+3jplzgq5Vg5myIBVw+Kyyp58IvdfLQxgcl9I/n0Htk4LYQ9CPByZfF9vZnSL4pPNicy8ZPtsl3OhkjA1VF6XinjP9jK2qOneeGW9rw4poMcBiCEHXE2Gph9S3vevK0z+1LyGDVvE5tOZOtdlqgF+SSug/0peYx59w9Sc0v4dEpP7ukXpXdJQogGMq57OD9PG4CfhwuTP93BGyuPySpLKycBd42WH8xgwodbcXUy8OMj/RjcOkjvkoQQDSy2mTfLpg3gtu7hzF8Xz10fbyczX1ZZWisJuKukqirvrovnkcV76BDmw9JH+9OqmbfeZQkhGom7i5HXx3fmrQmdOZSez6h5m1h3NEvvssRlSMBdBe20W/uZs/IYt3YNY/F9vQnwctW7LCGEDsZ2C2fZtAEEebty72c7mb3sMKYKueyONZGAq6WM/FImfLiVJXvTeHpEK96a0Bk3ZzlXnRCOrGWQF0sf7c+UflF8tuUUf3n3D45lFupdlqgiAVcLu5NyufmdPziZVcTHk3swbWisnHZLCAGAm7OR2be0Z+G9PckpKuOW+Zv5fMspVFUuvaM3CbgafLczhTs/2o6nq5Fxbof45+TRuLq6MmXKlGof8/nnn9O9e3eaNGlCeHg4zzzzDJWVlY1XtBAOav78+fTo0aPG9+jFhg4diqIodX6PDmkdxIrHr6Nvi6bMWnaY+z7fxZmisjotU9SNBFw1KswWZi87zDM/HqBXtD8/PdqfLm1ieO6555g6deoVH1tSUsLbb79NTk4O27dvZ82aNbzxxhuNVLkQjis0NLRW79FzFi9eXK9fPgO9XVk4pSezb27Hpvgcbnh7ExuOyzFzepHL5VxGWl4pj329l91JZ5naP5qZo9rgZDQwduxYAHbt2kVqamq1j3/44YfP/x4WFsbdd9/NunXrGrxuIRxdbd+jAPn5+bzwwgssWrSIvn371lsNiqIwpX80fVo05bGv93LPpzuY2j+af4xsLdeYa2QScP9j1eFM/v7DASrNFube0YUxXcLqvMyNGzfSvn37eqhOCFFfZs6cycMPP0xwcHCDLL9NcBOWTRvAK8uP8OkfiWxNOMPcO7rIYUWNSFZRVimrNDN72WEe+GI3Ef7u/PrYwHoJt4ULF7Jr1y6efvrpeqhSCFEfdu3axR9//MH06dMb9HncnI28MKYDn07pQVaBiZve2cxHG09itsgOKI1BAg5IzClm7Htb+GzLKe7tH8WPD/cjKsCzzstdunQpM2bMYMWKFQQEBNRDpUKIurJYLDzyyCPMnTsXJ6fGWYk1tE0zVj5xHUNaB/Ly8qPc8dFWks4UN8pzOzKHDjiLReXLbUncNG8TaXmlfDy5B7Nubl8v68l/++037r//fn7++Wc6duxYD9UKIepDQUEBu3bt4vbbbyc4OJiePXsCEB4ezqZNmxrseQO8XPlgYnfemtCZo5mFjJy7iU83J0o314AcdhtcQnYRM/57kB2JufRv2ZQ54zsT6ut+xcdUVlZSWVmJ2WzGbDZjMplwcnL607fAtWvXcvfdd7NkyRJ69erVkC9DCHGR2rxHfXx8SE9PP387JSWFXr16sXv3bgIDAxu0PkVRGNstnD4xTZm55CAv/hLHT/vSeGVsJ9qFNmnQ53ZIqqpeabI7ZRVmdf7aE2rss8vVjrN+U7/dmaxaLJZaPXbWrFkqcMk0a9YsNSkpSfX09FSTkpJUVVXVwYMHq0ajUfX09Dw/3XjjjQ35soRoPOtf0yYrVNv36MUSExNVQK2oqGjUWi0Wi/rTvjS1+0ur1Jh//qq+vDxOLTI1bg12otoMU9QrH21vN72zqqqsPHyaV1YcIelMCSM7BPPCmPYEebvpXZoQtmXD69rPQc/oW4edyCsp55XlR/l2VwpB3q78/YbWjOsWjsEgZ0uqpWoHyiEC7mBqPi/9GseOxFxig7x4dnRbubyNENdKAq5B7E46y0u/xLEvJY+OYT48N7otvWOa6l2WLXDMgNt5Kpf31sWz7lg2/p4uPHl9K+7oGSFX3BaiLiTgGozForJsfzqvrjhKZoGJgbEBPD4slh5R/nqXZs0cJ+DMFpX1x7L4YMNJdp46i5+HM/f2j2ZK/yiauDnrXZ4Qtk8CrsGVlFfyxdYkPtqYwJnicvq3bMr0obH0jvaXE73/mf0HXNKZYr7flcoPu1PJLDAR6uPG/dfFcHvPCDxcHHZnUSHqnwRcoykpr+Sr7cl8sCGBnKIy2oY04Z6+kYzpEoa7i5z2q4r9BZyqqsRnFbHmaBa/x51mV9JZDAoMahXIhB4RDGvbDBcnWRUpRL2TgGt0peVmlu5L4/MtpziaWUgTNyfGdQ9nXLdw2oc2cfSuzvYDTlVVUnJL2Z2cy65TZ9l0Iofk3BIA2oU0YVTHYMZ1DyfE58rHsgkh6kgCTjeqqrLz1Fk+33qKVYczqTCrtAzy4tauYdzUKYTIpnU/A5MNsp2AM1WYySooIy2vlIScIk5mFZOQU8Th9AKyC7VrK3m5OtEzyo9hbZsxtE1QjQdoCyHqkQScVcgrKWf5wUyW7k1jx6lcQLvC+LA2QQxtE0T3SD9H2aFO34DLyC/lww0JmC0qFlWbKswqJeWVFJoqKSqrpKC0gqzCMgpNl16byd3ZSEygJ62bedMt0o/ukX60auaNUY4REUIfEnBWJyW3hNVxp1l3LIttCWeoMKs0cXNiUOsgBsYG0CvKn8imHva6KlPfgDuaWcCED7ZiNCgYDQqKouBkUPB0dcLL1QlvN20K9HIlqIkbgd6uhPi40SLQi+AmbnLAoxDWRALOqhWVVbL5RDZrjmSx7lgWOUXlAAR5u9Iz2p/e0f70jPK3p0bBdlZR2qo333yT2bNnU1RUpHcpQjSo565zAeBfG8t1ruTqeXl5MXv2bJ566im9S2kUFovKyewitifmsvNULjsSc8nINwHa2rG2Id60D/WhQ1gT2of60KqZty3unCcB19BCQ0PJyMjQuwwhGpwtBxxASEjIJSdbdiSqqpJ6tpRdSbkcTC3gUHo+cekFFJVpm4acjQrN/T2IDvCiRaAnMYGexAR6ER3gSVNPF2tdxVltUXKAWD156qmnpIMTwsp5eXk5TPd2OYqiEOHvQYS/B7d21f5msagk5ZZwOD2fw+kFJGQXkZBdzMbj2ZSbLecf6+JkILiJmzb5uBHio21O8vVwwc/DGV8PZ3zcXfD1cMbL1QlXJ4PugSgdnBDi6sg2OIdQabZoe7NnF3PqTDGZ+SYyC0xk5Ju03/NNlwTg/1IUcHUy4OZsxM3JiJuzARcnA0aDgZs6hfDokJb1Vap0cEIIIWrPyWggsqlntcfWqapKfmkF+aUV5JVUkFdaQV5JOfmlFRSaKimrMGOqtFBabsZU9Xul2UKlRcXHvXFOmygBJ4S4Oj4RelcgrICiKPh6uODr4UKklV70QFZRCiGEsGXVrqK0uf1BhRBCiNqQgBNCCGGXJOCEEELYJQk4IYQQdkkCTgghhF2SgBNCCGGXJOCEEELYJQk4IYQQdkkCTgghhF2q9lRdL7zwgtPjjz/emLUIIYQQV2Xu3LlRQOqsWbMq//e+K52LMnzu3LkNVpQQQghRDxKBaODU/95xpYBLrXrQ1TyBaBgyvg1Pxrhhyfg2LEcf39TL/lVV1TpPs2fPVutjOTLJ+MoY2+ck4yvjq8dUXzuZvFBPyxGXJ+Pb8GSMG5aMb8OS8b2Mmi6XI4QQQtgkOUxACCGEXZKAE0IIYZdqFXCKoryhKEqioiiqoigdqplnhKIouxRFKVMU5Y36LdP+1XKM/09RlMOKouxXFGW3oig3NHadtqqW43uvoigHFEXZpyjKQUVRHmvsOm1Vbcb3onlbK4pSIp8TV6eW/4dnK4qSVfV/eJ+iKO82dp3WpLYd3FLgOiDpCvMkAPcDc+palIOqzRjvAHqqqtoZmAp8qyiKe2MUZwdqM74/Ap1VVe0C9AOeUhSlU2MUZwdqM74oimIEPqyaX1ydWo0xsEhV1S5V06ONUJfVutJxcOepqroZQFGUK80TXzXPmHqpzMHUcoxXXnTzAKAATanuGBBxXi3Ht+Cimx6AMyB7YdVCbca3ygzgF8CrahK1dBVjLKrINjjbNRk4qaqqhFs9UhTlFkVRDqN9S56jqupBvWuyF1Xd8A3Af/Suxc7dUbWqfZWiKH31LkZPEnA2SFGUQcBLwJ1612JvVFVdpqpqe6AVMElRlNZ612QPFEVxBj4GHlJV1ax3PXbsAyBaVdVOaJuLflIUpanONemmVqsohfWo+kb2JTBGVdVjetdjr1RVTVYUZQdwEyDjXHchQAtgedUqNl9AURSliaqqD+hamR1RVTXzot9XK4qSAnQANuhXlX4k4GyIoig9gW+B8aqq7tG7HnujKEobVVWPVv0eAAwB/qtvVfZBVdVkIODcbUVRZgNeqqo+rVtRdkhRlDBVVdOqfu8CROHAX9Bqe5jAPEVRUoFw4PeqbRQoirJcUZQeVb8PqJrnSeBBRVFSZTf22qvNGAPvAe7AhxftBtxRp5JtSi3H98GqwzD2AWuA+aqqrtKpZJtSy/EVdVDLMX5ZUZRDiqLsR1slPOnirs7RyKm6hBBC2CXZyUQIIYRdkoATQghhlyTghBBC2CUJOCGEEHZJAk4IIYRdkoATQghhlyTghBBC2CUJOCGEEHbp/wHE+zPfmRwxfQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_posterior(sensitivity_samples, var_names=['μ'], ref_val=np.log(4));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the original model for comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3wUdf7H8dfsbnolvZIQeq8BpIhUKQoIiBUL6lnRs5z6Uw/x9GyoZxcVxQYoeFKlSi/SIfQSSnoISUivuzu/Pya0E0iAJLPl83w89rFkd3bmM8Puvvf7ne/MKKqqIoQQQjgag94FCCGEEHVBAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk6IOqIoiqooSpPz/v5OUZQ39KxJCGciASeEEMIhScAJIYRwSBJwQgghHJIEnBB1y/28f/vrVoUQTkgCToi6db+iKEZFUToA/QEfRVFc9C5KCGcgASdE3fIEMoCvgYnAvUA/XSsSwkkocsFTIeqGoigq0FRV1US9axHCGUkLTgghhEOSgBNCCOGQpItSCCGEQ5IWnBBCCIckASeEEMIhmap5XvovhRBC2BKlphNKC04IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDqu5MJkI4DatV5WBmIZuP55BZUEZhmRlfdxfCfN1oGe5L2yg/PF3lIyOEvZBPq3B6OUXlzNySzPTNyWTklwHgajTg426ioKySSot2xjqTQaFrowAGtQrllk5R+Hm46Fm2EKIa1V0uR85FKRyWqqr8vDWFN38/QGG5md5NgxjRIZIejQOJ8PcAtFZddlE5e9Ly2XriNCsOnORIVhHebibu7NaQh6+PI9DbTec1EcKp1PhclBJwwikVllXy5MydrDp0im6NAnh9ZBuahfrU6LV70/L5et0xFiSk4+lq4pE+cTzYOw53F2MdVy2EQAJOiEtLyyvlge+2ciSriIk3tWJc9xgMhhp/Zs5KzCrk3SWHWLb/JLGBnrw+sg29mwbXQcVCiPNIwAlxMZn5ZYyZspH8kkq+uLszvZoGXfM81x/JZuK8vRzLLua+HrG8OKSFtOaEqDsScEL8r7ySCsZ++Sdpp0uZ+bfutIvyr7V5l1VaeHfJIb7dcJwWYT58dHtHmofVrMtTCHFFJOCEOJ/FqjLum81sO3Ga7+6Pp0eTa2+5XcyqQ1n8Y3YChWVmJt7ciru6xdTJcoRwYnLBUyHO99GKI2w8msMbt7Sps3AD6Ns8hMVPXU/3uEBenrOXf87di9lirbPlCSEuTQJOOLwNidl8svIIoztFMbZLdJ0vL9jHjW/vi+fh6+P4cVMS903bSn5pZZ0vVwhxIemiFA6tuNzMoP+sxc3FwMIJver9TCSztqbw8tw9RAd48u298cQGedXr8oVwQNJFKQTA+8sOk5ZXyjuj2+lymq2x8dH89EA3ThdXMGbKnxzKLKz3GoRwVhJwwmHtSslj2sbjjOseQ3xsgG51dIsLZPYjPTAocPtXf7I3LV+3WoRwJhJwwiGpqsqk+fsI9nbj+cHN9S6HJiHezHr4OjxdTdz59SZ2peTpXZIQDk8CTjikBbsz2JWSx3M3NsfH3TZOihwb5MUvD3fH39OVu6duZntSrt4lCeHQJOCEwymrtPDO4oO0CvdldKcovcu5QFQDT2Y9fB3BPm7cP20rh0/KPjkh6ooEnHA4P/6ZRFpeKa8Ma4nxKs4xWdfC/Nz5YXxXPFyN3PPNFtLySvUuSQiHJAEnHEpJhZkpa47Su2lQnR7Qfa2iAzz5fnxXiivMjPtmM7nFFXqXJITDkYATDuWnTUnkFFfwVP+mepdSrRZhvnxzbzxpp0sZ/91WSirMepckhEORgBMOo6TCzJdrjtG7aRBddDws4Ep0bRTAJ3d0ZHdqHv+YvZtqTrwghLgCEnDCYczckkJOcQV/H2D7rbfzDWodxotDWvD7ngy+XHtM73KEcBgScMIhmC1Wvl1/nK6xAXSOsY/W2/ke6h3HTe3CeXfJQdYePqV3OUI4BAk44RAW7c0kLa+Uh66P07uUq6IoCu+OaUezUB8mzNxJck6J3iUJYfck4ITdU1WVr9ceIy7Ii/4tQvQu56p5upr4clxnAP724zZKKyw6VySEfZOAE3Zvy/Fc9qTl80DvRhhs8Li3KxET6MXHd3Tk0MlC3vh9v97lCGHXJOCE3fthUxJ+Hi6M6mhbZy25Wn2aBfO36+OYvjmZZfsy9S5HCLslASfsWlZBGUv3ZnJr5yg8XI16l1Nrnh3YnDaRvrzw391kFZTpXY4QdkkCTti1GVuSMVtV7u4eo3cptcrVZODD2zpSWmnh2dkJWK1yfJwQV0oCTtitSouVGZuT6dMs2CGvlN0kxJuJN7Vm3ZFsvt1wXO9yhLA7EnDCbq04kEVWYbnNt95WrVpF37598fPzIzY29rLT7t+/ny5dutCgQQMaNGjAty+Pp4tfMe8uOcQRufKAEFdEAk7YrVnbUgjxcaNv82C9S7ksLy8vxo8fz+TJk6udNiIigl9//ZXc3Fyys7MZPnw4B6a/jqebkZfm7JGuSiGugAScsEvRMTHM/eELUqY+jp+vDw888AAnT55kyJAh+Pj4MGDAAE6fPg3Apk2b6NGjB/7+/rRv357Vq1efnc+0adNo2bIlPj4+xMXF8eWXX559bvXq1URFRfH+++8TEhJCeHg406ZNu+Jau3btyrhx44iLq/4gdH9/f2JjY1EUBVVVMRqNHD92lJeHtmTridP8vDXlipcvhLOSgBN2qaTcQvHBjcz/fTGHDx9mwYIFDBkyhDfffJPs7GysVisff/wxaWlpDBs2jFdeeYXc3Fzee+89Ro8ezalT2umwQkJCWLhwIQUFBUybNo2nn36aHTt2nF1OZmYm+fn5pKWl8c033/D444+fDc63334bf3//S96uhb+/P+7u7kyYMIGXXnqJMZ2juC4ukLcWH5BRlULUkAScsDuqqlJcYabD4NuJbxVHZGQkvXv3plu3bnTs2BE3NzduueUWdu7cyU8//cTQoUMZOnQoBoOBgQMH0qVLFxYtWgTAsGHDaNy4MYqi0KdPHwYNGsS6devOLsvFxYWJEyfi4uLC0KFD8fb25tChQwC8+OKL5OXlXfJ2LfLy8sjPz+fTTz+lY8eOKIrCm6PaUm62MmnBvmuatxDOQgJO2J3Nx3MxW1QGd2l+9jEPDw9CQ0Mv+LuoqIikpCRmz559Qctq/fr1ZGRkALB48WK6d+9OQEAA/v7+LFq0iOzs7LPzCQwMxGQynf3b09OToqKielhLbd/dI488wj333ENWVhaNgrx4qn9TFu3JZPn+k/VSgxD2TAJO2J1ZW1MwKNq11KoTHR3NuHHjLmhZFRcX8+KLL1JeXs7o0aN57rnnOHnyJHl5eQwdOrTG12R788038fb2vuStNlitVkpKSkhLSwPgb9fH0TzUh4nz9soFUoWohgScsCsFZZUs2puBp6sJN5fqz1xy9913s2DBApYuXYrFYqGsrIzVq1eTmppKRUUF5eXlBAcHYzKZWLx4McuWLatxLS+99BJFRUWXvJ1htVopKyujsrISVVUpKyujoqLiovNcvnw5O3fuxGKxUFBQwDPPPEODBg1o2bIlAC5GA2/c0oaM/DKmrJFrxwlxORJwwq4sSEinrNKKp1vNTssVHR3NvHnzePPNNwkODiY6OprJkydjtVrx8fHh448/ZuzYsTRo0IAZM2YwfPjwWq957dq1eHh4MHToUJKTk/Hw8GDQoEFnn2/dujXTp08HtH1vd9xxB35+fjRu3JjExESWLFmCu7v72enjYwO4uX0EX645Slpeaa3XK4SjUKrpjpGDboRNGfHpesrNVhY/1RtFse8rB1yLtLxS+r+/mgEtQ/n0zk56lyNEfarxB19acMJuJGYVkpCaz5jOUU4dbgCR/h481DuOhbsz2JuWr3c5QtgkCThhN+buTMegwPAOEXqXYhMeuj4OPw8X3l92SO9ShLBJEnDCLqiqyryENHo2CSLEx736FzgBX3cXHu4Tx6pDp9ielKt3OULYHAk4YRd2JJ8mJbeUkR0i9S7FptzXI5YgbzcmLz1U48MbhHAWEnDCLszdmY67i4Eb24TpXYpN8XQ18Xjfxmw6lsuGxBy9yxHCpkjACZtXabHy+54MBrQMxdvNVP0LnMyd3RoS4efO5GXSihPifBJwwuatO3KK3OIK6Z68BDeTkSf7NyUhJY8/DmTpXY4QNkMCTti8uTvT8fd04fpmtn3dNz2N6RxFbKAnH/5xWFpxQlSRgBM2rbjczPL9JxnWNhxXk7xdL8VkNPDoDY3Zl14g++KEqCLfGMKmLdufSWmlhZEdpXuyOiM7RhLi48aUNUf1LkUImyABJ2za3J3pRPp70LlhA71LsXluJiPjezVifWI2e1Ll7CZCSMAJm5VdVM76xGxGdIjAYHDuU3PV1J3dGuLjZmLKWmnFCSEBJ2zWwoR0LFZVuievgK+7C3d1j2HxngyScoovP/GumdpNCAclASds1txd6bQM96VZqI/epdiV8T1jMRkMfL2umuvF5adoNyEclAScsEknsovZlZLHSDmx8hUL8XVnVKdIZm9LJbuoXO9yhNCNBJywSfN2paPIlQOu2kPXx1FhsfL9xhN6lyKEbiTghM1RVZV5u9Lo1iiAcD8PvcuxS42DvenfIoSZW5IpN1v0LkcIXUjACZuzJy2fY9nFcmqua3Rvj1iyiyr4fXeG3qUIoQsJOGFz5u5Mx9VoYEjbcL1LsWu9mgTRONhLuimF05KAEzbFYlVZsDudvi2C8fNw0bscu6YoCvf2iCUhNZ+dyaf1LkeIeicBJ2zKxqPZnCosl+7JWjKqUxTebiZpxQmnJAEnbMrcnen4uJvo2yJE71IcgrebiTGdo/h9TwZZhWV6lyNEvZKAEzajrNLC0n2ZDGkThruLUe9yHMY918VQaVH5eYsc1C2ciwScsBl/HDhJUblZuidrWVywN9c3C2b65iQqLVa9yxGi3kjACZsxd2c6ob5udIsL1LsUh3NfjxhOFpSzdF+m3qUIUW8k4IRNyCupYM3hLIa3j8AoVw6odTc0CyHS30O6KYVTkYATNmHRnkwqLSojpHuyThgMCrfFR7M+MZvknBK9yxGiXkjACZswd1caTUK8aR3hq3cpDmtM5ygMCszaJq044Rwk4ITuUnJL2HI8l1s6RqIo0j1ZVyL8PejTLJjZ21Mwy2AT4QQk4ITu5iekAzC8vVw5oK7dFt+QkwXlrDl8Su9ShKhzEnBCV6qq8tuOVLrGBhAd4Kl3OQ6vf8sQgrzd+HmrdFMKxycBJ3S1L72Ao6eKGdlRBpfUBxejgTGdo1h5MIuicrPe5QhRpyTghK7m7EzD1WhgmFw5oN7cFh+NxaqyLz1f71KEqFMScEI3ZouV+QlVVw7wlCsH1JdGQV50axTA3rR8VFXVuxwh6owEnNDNxqM5nCos5xbpnqx3t3eNJq+kkpTTpXqXIkSdkYATupm7Mw1fdxM3NJcrB9S3wa3DcTUZ2J9eoHcpQtQZCTihi5IKM0v2ZTKsXbhcOUAHHq5Gmof6cORkISUVMthEOCaT3gWIepSxG/54FdJ3QmUZ+DeErg9pN4C9v8HqtyE/FVC15+MfPPf8xZjLYflEOLAQirPAIwAa94Mb/w2eAZB7HOY8Apm7Ibw93PIlNIhh2b6TPGz9hQfSDoN1HRhq4bfWH5PgwALISdT+HvE5dLzr8q/5YQRk7oWyfPDwh4bdYdAb0CBWe74kF5a+DEdXQmkueIVAy5tg4L/A5KZNU1kKq/4Ne+dA0UnwDISOd0P/f2rPz58AyZu17WpyhcguMOh1CGlZ9foyWPsu7JkNhSchIA76/h+0GlGn26RlhC970vL588+19E/9ApL+BKtZW/eb/gMx1118npumwO5fIPcYWCohqCn0eQFaDNWeT98J8yZoz8f2glumaO8F0N4L5YVw+/RrXzchqiEtOGfy813aF7VPBDQdCNmHYdFzcHyt9nx+CvhHQ/vbILY3nDp44fMXs+4D2DxF+9JqNQJQIWEGLH1Je375P7VwazYY0ndpfwPrt27jYdNCPEd8ePXhVpILqdvO/Z26DfyitZCtqYJ0aDJA+9I3mLQwmPvYueeXvqytDyq0Gqmt5+YpsO597XlVhV/uho2fgNEEHe7QQjL32Ll57PgB3Hyg7WjtPnE5/DhKCzbQttW698HgAu1vh6JMmHUvpGypQf0ZkLnn0s9fZptE+nsQ7VFKj9V3wZFlENVZW75XEBSkXXqeBxZA6WloPhRCW0HGLpg1DjIStOcXPAUFqdp7LHE5rP9Aezx5E+yfB4Pfqn69hKgF0oJzFpZK7UsHYPRU7Yvpyz7al1NesvZ4z6e02xmf94CsfXA6CRpdYr6nj2v3ncZprbbNX8Hif5yb56nD0Oh6uHUaTL8Vsg5yqrCcISn/4XDYYNo1jL+y9bBaIHEF7JoOhxZDm9EQ1UV77r6F2v0XvbTWVk08sfXcv/fP176oTyf9df16Pwfd/qYF3p+fnlu/42sg8Q8IagYPrwMX978uY/wyaNitan5J8FE7KEzXfkBEdID9c7Xnhn+stXiCW8DS/9NC785f/jo/czkc/B12zdB+sFz/Dwhre/H1u8w2MSgKN7gexqOkhMLuz+IzeGI1G6vKgEkQ2Vn7YWK1wCed4PQJOL5Oa6WfOqy9H4ZOhs+6QdZBbbrfn4Nez2g9A0LUAwk4Z2F0gW6PwqbP4LeHtG6wjAQIbQstbjo3Xep22DNL61rM2gdBzaHFsEvPt8t4OLgIdvwIRVnaF76LJ/R4Uns+uBkkroRfH4AT66HpQHb98TPxhsPkDPu+5vWfOqyFWsLPWgsnpBXc8CK0u+3qtsf5Vr+jBc7hZaAYoeeT557r9ojWClr3HqRuhcNLtdbQmW7bY2u0e1cvmNITCjMhvAMMeftc6JwJN9B+aAAoBvAJ0/5tqgrFjAQtOE7u1f7O3HthnWk7tG2w51coy9OmHfwWtL31qlc9qFI7TVru4c34JMRqrcjWI2HAa+B6iTPLRP/Pj5Iz6+Rbdaq14Gaw979QnK31EjQdCFunQmXxhdtWiDomAedMWgyDgwu0L9CTe7UvsxbDtG6zM04d1LrgQPsSbjIAXL0vPc/g5tC4LxyYrwUjaN2bZ/YvDXwdik7BoUXar/u+L9Pmi+FM97qXx0+vh0V3aNNd95i23+pipo+FI0vBJxzajtFCLbzdtW2L8+38CfKrWmRBzbU6z4jsDNHdIGn9ufVrNQIaVDVpS3K0+/Sd0HyYVuOJdVrNE7ZpwXdGeRHMfbRqfZ84F3C9n4Hfn9W6Ks907YK2P++Mz7rDqQPacrs9rG2DwMbXvOoGs3aYQEDuTtSOo1ES/4AtX2lBP+Tt6mew9CWtOzO6G7Qcrj1280faPrjDS6HJQOh4D3wzAEZ/A39+Bgkzte1yw0vQbNA1r4MQlyL74JxFSS5MH6N1rd2/BF44obUw1rwN2749N13Hu2DiaZiwQ3t+02ew8aNLz3fh01q4xT8IL2dq3Vcn1sHs+7TnAxrBA0vh5QwYv4TczT+TbXYntG0/mP8EdH0Q4sdrAzFOHb74MvKrulaDm0NoG22etenpPfBSOgx7H7IPwYyxUFGsPTf7Xi3c+r+qrV/8g9p+pIVPa897BVXV1gLumAHj5oKbr9YiTN91bhnFOfD9zZC6BTrdqw1SOSP+QXhwJfR9Bfq+DMM/uXDeZ7eBAqGttW3gF1U76+6itdI+qRzO7k6vw+A3tccPLbr866wWmPeE9mMooqPWlWqs+r0c0REeXQ8vp8Nds2D9fyCml/aDacVrMOjf2kCk2fdp+/KEqCMScM7i9AmoLNFabZGdwKOBFhigdSOBNoACtH0rgY21LyqAnKPn5nPqsHYzl2t/Zx3U7iM6gouH1uI5f57/U4P3js951Xw/A4NyQbVq00d21v6dte/itT+8Bm77CUweWihObgqz79f2wZ3pHquJ4hyt9oL0qvUt0gaJgNaiONNVW5YPxacuXL+oLtr6ndkmZ9YvtPWll3em9ZaXDN8OgvQd0OtpbV/b+ZcFMldoAzz6/AP6PK+NZgSIu+HcNM8d1kYjluXDrHu0bTDvca2L1HoNl77x1o5BNBoM/LYj9bztUdVqt1Se+z8/s5zKMvhlHOz8EeL6wr0LtffTxSRv1vYxDn7r3GCY6HitlVxZfOFgHCFqmXRROovg5tqXUOlp+H641gra86v2XMPu2v2XfaBBjDZMvCBD6xYEaNz/3Hw+q9r/8vA6rZuwYTet6+yPSdo+qmOrL5znedTFL7CEnvg27YFftIf24O/PAVVfqoFNL1670QVa3qzdirK0Ieq7ZsDM26HjOBjxqTbdug8g+4g2GhS00Ysn1kOne7Qh71u+0lqszYdpra3987Th+VFdtfA6tupcHf4xVevRTRvIMX+CFjgHf79w/VpWdVeeOggz74SKQigv0Fq/oW20ab4ZBIUZ2mjGyjJY/KL2eNtbtWDb8T3snqUN/Mk6ACmbwc1PGzxyhqunNsKx/e3aQJWEmdpt509ay6/PedOe73LbBKDhdZB7jCdc5rN0ZyZq4gEU0EaDgvZj4Mz/+QtJ2qEU85+AQ79r+w4DG8PKN7TnIztDu/P2B1otsOhZLdQbxGgDcQB+Ha/N1+gK/rEXr1uIWiAB5yxcveCuX2Hl69pghowEbaBJl/u1kYigfYEfWa59Abp4aq2V+Acv/NL6X4Pe0FqFh5dooeMRAO3vuLALDuDwUszHNzKp5F1e7xwNYeFad+bGqnAaMAnC2lS/Ht4h0GOCdkvbceFw9sQVWnfiGSmbtFtsr4sf0xXYBLyCtSCvLNPm3fFu6PPiuRbWyC9g+atacO+aoR0HF/+QVi9o3XJ3/xcWv6AFpIuntj0HvXGuy64wQ7vPT4HNX5xbflhbLeACGmk/PHbN1L70mw3W5n+pfWwNYrQBNn1e0P6vrJc5UPty2wS0HzO3fAV/vM2wgrWUqJF4DXwduj9+6XkWVK2PuUwbPHJG+zsvfK9s+1ZrJZ8Zmdt8iDbQKWEGuHjBzR+DV+CllyPENVKqOdmqnIlV1JonZuxgfWI2m1/qj5tJzl6iuzXvavd9nsdssXLd2yvpGO3PV/d00bcuIS5PqX4SjeyDE/Uiv6SSZftPMrJDpISbDTIZDYzsEMGqQ1nkFlfoXY4QtUICTtSL+QlpVJitjOlcS6P/RK0b3TmKSovKgoR0vUsRolZIwIl6MXt7Ki3DfWkT6ad3KeISWoT50ircVxtNKYQDkIATde5QZiG7U/O5VVpvNm905ygSUvM5crJQ71KEuGYScKLOzd6WgotRYaRc2NTmDW8fgdGg8NvOy5xsWQg7IQEn6lSlxcrcXWkMaBlKgJer3uWIagT7uNGnWTBzd6ZhscogamHfJOBEnVp5MIvsogpu7SLdk/ZiVKdIMvLL2HQsR+9ShLgmEnCiTs3ckkyYrzvXNw3WuxRRQwNahuLjbuK/MthE2DkJOFFn0vJKWXP4FGPjozEZ5a1mL9xdjNzULpwlezMpLr/MWVKEsHHyrSPqzC9btfMf3hYfrXMl4kqN6hRFSYWFpfsy9S5FiKsmASfqhNliZdbWFPo0CybS30PvcsQV6hLTgOgAD37bIaMphf2SgBN1YvWhU2QWlHFH14Z6lyKugqIojOoYxYaj2WTkl+pdjhBXRQJO1ImftyYT4uNGvxYhepcirtKoTpGoKszdKafuEvZJAk7Uuoz8UlYezOLWLlG4yOASuxUT6EWXmAb8tiOVaq46IoRNkm8fUetmbU3FqsLt8dI9ae9Gd47iSFYRe9MK9C5FiCsmASdqlcWq8svWZHo3DSI6wFPvcsQ1Gto2HFeTQY6JE3ZJAk7UqtWHskjPL+NOGVziEPw8XBjYKpT5CelUWqx6lyPEFZGAE7Xqx01JhPq6MaBVqN6liFoyulMkucUVrDl0Su9ShLgiEnCi1iTlFLPm8Cnu6NpQBpc4kN5NgwnyduW3ndJNKeyLfAuJWjN9czJGRZFj3xyMi9HA8PaR/LE/i/ySSr3LEaLGJOBErSirtDBrWwo3tg4j1Ndd73JELRvVKZIKi5WFe+SYOGE/JOBErViQkE5eSSV3d4/RuxRRB1pH+NI81EdO3SXsigScqBU/bUqiaYg33eMC9C5F1AFFURjVKZLtSac5nl2sdzlC1IgEnLhmCSl5JKTmM+66GBRF0bscUUdGdozEoMAcOSZO2AkJOHHNftyUhJerkVs6RupdiqhDob7u9GwSxH93pGG1yqm7hO2TgBPX5HRxBQsS0rmlUyQ+7i56lyPq2Ngu0aTllbLxaI7epQhRLQk4cU1mb0+h3GyVwSVOYmCrUPw8XJi1LUXvUoSolgScuGpWq8pPm5LpGhtAizBfvcsR9cDdxcjIDhEs2Zcpx8QJmycBJ67a2iOnSM4tYdx10npzJrd2iabCbGV+ghwyIGybBJy4aj/+mUSQtxs3tg7TuxRRj9pE+tEq3JdZ22Q0pbBtEnDiqiTlFLPyUBZ3dI3G1SRvI2cztksUe9Ly2Z8u14kTtku+mcRV+W7jCUwGRQaXOKmRHSNxNRqYvV0GmwjbJQEnrlhhWSWzt6VyU7sIOe+kk/L3dGVQ61Dm7kyj3GzRuxwhLkoCTlyx2dtSKSo3c3/PWL1LEToa2yWa0yWV/LE/S+9ShLgoCThxRSxWle82nqBLTAPaRfnrXY7QUc8mQUT4ufPz1mS9SxHioiTgxBVZeTCL5NwS7u/ZSO9ShM6MBoXb4huyPjGblNwSvcsR4i8k4MQV+Xb9cSL83LmxdajepQgbMDY+CgX4ZasMNhG2RwJO1NiBjAL+PJbDPT1iMRnlrSMg3M+Dfi1C+GVbCpUWq97lCHEB+ZYSNTZtw3E8XIzcHh+tdynChtzRtSGnCstZcUAGmwjbIgEnaiSnqJy5u9IZ3TkSf09XvcsRNqRPs2DC/dyZuUUGmwjbIgEnamTG5mQqzFbu6yGDS8SFTEYDY7tEs/bIKXxIHLEAABy3SURBVBlsImyKBJyoVoXZyg+bkujTLJgmId56lyNs0Nj4aBSQy+gImyIBJ6r1+550ThWWy4Hd4pIi/T24oXkIv2xNwSyDTYSNkIATl6WqKlPXHadxsBfXNw3Wuxxhw+7o2pCswnJWHJTBJsI2SMCJy1qfmM2+9AL+dn0cBoOidznChvVtrg02+WlTkt6lCAFIwIlqTFlzlFBfN0Z2jNS7FGHjTEYDd3VryLoj2Rw9VaR3OUJIwIlL25Oaz4bEHMb3bISbyah3OcIO3BbfEBejwo9/SitO6E8CTlzSlDVH8XE3cWe3hnqXIuxEsI8bw9qG89/t2hUnhNCTBJy4qBPZxSzem8Hd3WPwcXfRuxxhR+7pEUthuZk5O9P0LkU4OQk4cVFfrTuGyWiQQwPEFesY7U/bSD9+2HgCVVX1Lkc4MQk48RdZhWX8uj2V0Z2iCPGRK3aLK6MoCvdcF8ORrCL+PJajdznCiUnAib+YtuEElRYrf7s+Tu9ShJ26uX0EDTxdZLCJ0JUEnLhAYVklP21KYkibMBoFeeldjrBT7i5GxsZHs2z/SdLzSvUuRzgpCThxgRmbkyksM/NIn8Z6lyLs3N3dYrCqKjM2y1UGhD4k4MRZ5WYL36w/Ts8mgbSL8te7HGHnogM86d8ilJlbkik3W/QuRzghCThx1tydaWQVlkvrTdSae3vEkFNcwaI9GXqXIpyQBJwAoNJi5bNVR2kb6UevJkF6lyMcRM/GQcQFe/H9RhlsIuqfBJwAYN6udJJzS3iqf1MURU6qLGqHwaBwT/cYdqXksSslT+9yhJORgBOYLVY+XXmEVuG+9G8Zonc5wsGM7hyFj5uJqeuO6V2KcDIScIIFu9M5kVPCk9J6E3XAx92FO7s1ZPHeTFJyS/QuRzgRCTgnZ7GqfLIykRZhPgxqFap3OcJB3dczFgXtJAJC1BcJOCe3cHc6x04V82T/pnJBU1Fnwv08uLl9BL9sTSa/tFLvcoSTkIBzYtaq1luzUG8Gtw7Tu5waOXDgAP369cPPz48mTZowZ86ci0732muvoSgKf/zxx9nHZsyYQXh4OI0aNWL16tVnHz969Cg9evTAYrn0sVrfffcdvXr1+svjsbGxZ5fx3XffYTQa8fb2xtfXlw4dOrBw4UIAVq9ejcFgwNvbG29vb6Kiohg7dixbt269ms1glx7s3YjiCgszt8iB36J+SMA5sYV7MkjMKmJCP/tovZnNZkaMGMFNN91Ebm4uX331FXfffTeHDx++YLqjR4/y66+/Eh4efsFrX3zxRXbs2MEnn3zCE088cfa5J598kg8++ACj8dov6nrddddRVFREXl4eDzzwAGPHjiU3NxeAiIgIioqKKCwsZNOmTbRo0YLevXuzYsWKa16uPWgd4UfPJoFM23CcCrNV73KEE5CAc1Jmi5UPlx+meagPQ9uGV/8CG3Dw4EHS09N5+umnMRqN9OvXj549e/Ljjz9eMN0TTzzBO++8g6ur69nHcnJyiIyMJDw8nAEDBnDsmDai79dffyUyMpLu3bvXaq0Gg4Hx48dTWlp6dllnKIpCVFQU//rXv3jwwQd54YUXanXZtuzB3nGcLChn4e50vUsRTkACzkn9tiONY9nFPDOoGUY7aL0BF722mKqq7N279+zfs2fPxtXVlaFDh14wXXBwMDk5OaSmprJ8+XJat25NUVERb7zxBm+99Vat12o2m5k6dSre3t40bdr0ktONGjWKHTt2UFxcXOs12KIbmgXTNMSbr9Yek2vFiTonAeeEys0WPvzjMO2j/Oxq5GSLFi0ICQlh8uTJVFZWsmzZMtasWUNJiTb0vKioiJdeeokPP/zwL681GAx88cUXjBkzhvfee4+vv/6aiRMnMmHCBPbs2UPfvn258cYbLwjL/7Vp0yb8/f0vuCUnJ190mrCwMGbOnMmcOXPw8/O75DwjIiJQVZW8POc4CFpRFB66Po6DmYWsT8zWuxzh4Ex6FyDq34zNyaTnl/HumPZ2ddybi4sLc+fOZcKECbzzzjt06dKFsWPH4ubmBsCrr77KuHHjaNSo0UVf379/f/r37w/A7t272bZtG5MnTyY2Npb169eTkpLCgw8+yKZNmy76+u7du7N+/foLHouNja12mstJS0tDURT8/Z3n5NYjOkQweekhvl53nN5Ng/UuRzgwacE5mZIKM5+tSqR7XAA9mwTqXc4Va9euHWvWrCEnJ4elS5dy7NgxunbtCsCKFSv4+OOPCQsLIywsjJSUFMaOHcs777xzwTxUVeWJJ57g448/Jjs7G4vFQkxMDPHx8ezevbte12fOnDl06tQJLy/nufaem8nIfT1iWXv4FAczC/QuRzgwacE5mWkbTpBdVMGX45rbVevtjN27d9OsWTOsViuff/45GRkZ3HfffYAWcJWV546xio+P54MPPmDIkCEXzGPq1Kl07NiRDh06YDabKS0tZf/+/SQnJxMXV/dXMVdVlfT0dKZOncrUqVOZP39+nS/T1tzVrSGfrkxk6rrjvHdre73LEQ5KAs6J5JdU8uWao/RvEULnmAC9y7kqP/74I1OnTqWyspLevXuzfPnys12UgYEXtkiNRiMNGjTA29v77GPZ2dl89NFHbNy4EQCTycSnn35Kv379cHd3Z9q0aXVWe3p6Ot7e3qiqip+fHz169GD16tW1PoLTHvh7ujK2SxQztiTz7KBmhPt56F2ScEBKNSOZZJiTA3lj4X6+2XCc3yf0plWEr97lCL2teVe77/O8LotPyS3hhvdWc+91sUy8uZUuNQi7VOOuJ9kH5ySScor5/s8TjO0cLeEmbEJ0gCcjOkQwc0syucUVepcjHJAEnJN4e/FBXIwGnh3UTO9ShDjrsRsaU2a28N2G43qXIhyQBJwT2HI8l8V7M3mkT2NCfN31LkeIs5qE+HBjqzC+23iCwjI5CbOoXRJwDs5qVXnj9/2E+brzUO+6HyEoxJV6rG9jCsrMTN8sJ2EWtUsCzsHNS0hjd2o+zw9ujofrtZ9MWIja1i7Kn95Ng5i67jhllZe+ooMQV0oCzoGVVJh5d8kh2kb6MbJDpN7lCHFJj93QhOyicmZvS9G7FOFAJOAc2CcrE8nIL2Piza3s4nI4wnl1jwugU0N/pqw5JpfSEbVGAs5BJWYVMXXdMUZ3iiI+1j4P6hbOQ1EUnuzflLS8UmZvl1acqB0ScA5IVVVenb8XDxcj/ze0hd7lCFEjfZoF07GhP5+tTKTcLPvixLWTgHNAC3dnsCExh3/c2Jwgbze9yxGiRhRF4ekBzUjPL2PWVmnFiWsnAedgisrNvPH7ftpE+nJntxi9yxHiivRuGkSXmAZ8tuqojKgU10wCzsF8uPwwWYXlvD6ijd1cqVuIMxRF4emBzcgsKOMXacWJayQB50ASUvL4dsNxbo9vSMeGDfQuR4ir0qNxIF1jA/hsVaK04sQ1kYBzEOVmC8/NTiDU110Glgi7pigKfx/YlKzCcmbI2U3ENZCAcxAfrzjCkawi3hzVFl93F73LEeKa9GgcxHVxgXy2KpGicrPe5Qg7JQHnAPak5jNlzTFu7RxF3+YhepcjRK14YUgLcoor+GrNUb1LEXZKAs7OVZitPDc7gSBvV165SS4aKRxHh2h/hrUL5+t1x8kqKNO7HGGHJODs3Ccrj3DoZCFvjWqLn4d0TQrH8vyNzTFbrfznjyN6lyLskAScHdt6IpfPViUypnMU/VqE6l2OELUuJtCLu7rFMGtbColZRXqXI+yMBJydyi+t5O8/7yKqgSeThrfWuxwh6syEfk3wcDHy7pKDepci7IwEnB1SVZV/zt1LZkEZH93eAW83k94lCVFnAr3deKRPHMv2n2TriVy9yxF2RALODs3Zmcb8hHSeHtBUDugWTuGBXnGE+brz2oJ9WKyq3uUIOyEBZ2eSc0qYOG8fXRsF8OgNTfQuR4h64eFq5OVhLdmbVsCMLXLwt6gZCTg7UlZp4fEZO1AU+M9tHeRck8Kp3NQunB6NA5m85CA5ReV6lyPsgAScHfnXwv3sScvn/VvbE+nvoXc5QtQrRVH414jWlFRYeEcGnIgakICzE//dnsqMzck80qcxg1qH6V2OELpoEuLDA70bMWtbKtuTZMCJuDwJODtwMLOAl+fuoXtcAM8NaqZ3OULo6sl+TQn3c+efc/dhtlj1LkfYMAk4G1dQVsmjP+3A192Fj+/oiMko/2XCuXm5mfjnTa3Yn1HAdxtP6F2OsGHybWnDLFaVv/+8i+TcEj69sxMhPu56lySETRjSJowBLUN4d+khErMK9S5H2CgJOBv2zpKDrDyYxaThrenaKEDvcoSwGYqi8Oaotni5GnlmVgKV0lUpLkICzkbN3pbCV2uPMa57DOO6x+hdjhA2J8THnX/f0pbdqfl8vkouqSP+SgLOBm07kcvLc/bSs0kgE2+WS+AIcSlD24YzokMEn6w8wp7UfL3LETZGAs7GpOSW8PCP24ls4MHnd3bGRQaVCHFZ/xrehkBvV56ZtYuySove5QgbIt+eNiSnqJx7vt2C2aoy9d4u+HnK9d2EqI6fpwvvjmnPkawiJs7bi6rKuSqFRgLORpRUmBn//TbS80r55t4uNA721rskIexGn2bBPNmvCbO2pfLTZjlXpdBIwNmASouVx6fvYE9qHp/c0ZEusTJiUogr9fcBzejbPJjX5u+Ty+oIQAJOd6qq8tJve1h16BRvjGwrp+ES4ioZDAof3t6RqAYePPrTDjLzy/QuSehMAk5Hqqryr4X7mb09laf6N+XObg31LkkIu+bn4cJX93ShpMLMIz9tl0EnTk4CTieqqvLOkkNM23CC8T0b8fcBTfUuSQiH0CzUhw/GtichNY/Hpu+gwiwHgTsrCTidfLTiCFPWHOWubg35500tURS5tpsQtWVwm3DeGNmGlQezeGbWLrkKuJMy6V2AM5qy5igf/nGEMZ2jeH1EGwk3IerAXd1iKCwz8/big3i7mXhrVFv5rDkZCbh6pKoqn6xM5IPlh7m5fQTvjG6HQa7KLUSdeaRPYwrLKvls1VG83Ey8Mkx6S5yJBFw9UVWVtxYf5Ku1xxjVKZJ3R7fDKOEmRJ17blBzisstfLP+OEVlZv59Sxu57JSTkICrB1aryivz9jJjczLjusfw2vDW0nITop4oisKrN7fCx93EJysTySku55M7OuHhatS7NFHH5GdMHSs3W3h61i5mbE7m0Rsa868REm5C1DdFUXh2UHP+NaI1Kw5mcdfUTeSVVOhdlqhjEnB1KK+kgnHfbGHernSeH9ycFwa34LPPPqNLly64ublx3333XfK133//PZ07d8bX15eoqCief/55zGZz/RUvRD369NNPa/S5OF+/fv1QFOWKPhf3XBfLZ3d2Ym9aAaM+3ygXS3VwEnB1JCmnmFFfbGRXch4f3d6Bx25oAkBERASvvPIK48ePv+zrS0pK+PDDD8nOzmbz5s2sWLGC9957rz5KF6Le1fRzccb06dOv+gff0Lbh/PRgNwrKKhnx6QaW7M24qvkI2ycBVwe2J+Vyy+cbyS2u4KcHuzGiQ+TZ50aNGsXIkSMJDAy87DweffRRevfujaurK5GRkdx1111s2LChrksXQhc1/VwA5Ofn89prr/Huu+9e9fK6NgpgwYReNA314ZGfdvDOkoNyrJwDkoCrRaqq8sOfJ7j9q034uJv47dEedG1UOydOXrt2La1bt66VeQlhz1566SUeffRRwsKu7byt4X4e/PJwd+7o2pAvVh/l7qmbSc8rraUqhS2QgKslJRVmnv5lFxPn7aNXkyDmPd6TuFq65M20adPYtm0bzz33XK3MTwh7tW3bNjZs2MCECRNqZX5uJiNvjWrLu2PakZCax+AP17IgIb1W5i30JwFXC46eKuKWzzYyLyGdZwY245t74/H3dK2Vec+dO5cXX3yRxYsXExQUVCvzFMIeWa1WHnvsMT766CNMpto9wmlsl2gWPdmbuGBvJszcydO/7KKgrLJWlyHqnwTcNbBaVaZtOM6wj9eRVVjGd/d35cn+TWvtMIAlS5bw0EMPsWDBAtq2bVsr8xTCXhUUFLBt2zZuu+02wsLCiI+PByAqKop169Zd8/xjg7z49ZHr+PuApsxPSGfgB2tYvv/kNc9X6EcO9L5KqadL+Mfs3fx5LIe+zYN5e3Q7Qn3dq32d2WzGbDZjsViwWCyUlZVhMpn+8ot05cqV3HXXXcyZM4euXbvW1WoIYRNq8rnw8/MjPf1c92FKSgpdu3Zl+/btBAcH10odJqOh6sKpIbzw39089MM2hrULZ9LNrQn2cauVZYh6pKrq5W7if1SaLep3G46rrScuUVv9c7E6c3OSarVaa/z6V199VQUuuL366qtqUlKS6uXlpSYlJamqqqo33HCDajQaVS8vr7O3wYMH19VqCWe0+h3tZgNq+rk43/Hjx1VAraysrJOaKswW9ZMVh9WmLy1S201aqk7flKSaLTX/rIs6U11unb0pqnrZobEybvY8m47lMGn+Pg5mFtKrSRBvjWpLdICn3mUJcXXWVA2z7/O8vnXYuMSsIl6es4fNx3NpHeHLpOGtiY+tndHR4qrUeB+QBFwNJOeUMHnZIRYkpBPp78E/b2rJja3D5Kzkwr5JwNWYqqr8vieDf/9+gIz8Moa3j+D/hrYg3M9D79KckQRcbUjJLeHTlYn8uiMVk0HhkT6NeaRPYzlJq3AMEnBXrKTCzJTVR5my9hhGReHxvo0Z36sRnq4ynKEeScBdi6Onipi67jizt6VgMCjc2bUhj93QmJAaDCIRwm5IwF21lNwS/v37AZbsyyTI243HbmjMnd0a4u4iP37rgQTclVJVlfWJ2Xy7/jirDp3C1WjgtvhoHuvbWLohhGOSgLtmW0/k8sGyw/x5LIcwX3ce69uYWztHSy9P3ZKAq6msgjLm7Exj9vZUErOKCPJ2Y1z3GO7s1lCGBQvHJgFXazYezeaDZYfZlnQaf08X7uzakHuuiyXMT3p96oAE3OUUlFWy6mAWc3amsfbwKawqdI5pwO3x0QzvEIGbSX59CScgAVerVFVl64nTfLP+GMv2n8SoKAxuE8aYzlH0ahIkVxGvPRJw/ystr5SVB7NYti+TTcdyqLSohPu5M6pTJKM7RdXaeSOFsBsScHUmOaeE7zae4LedqeSVVBLi48bIjpEMbx9B6whfGYF9bZw74FRVJT2/jB1Jp9l4NIeNR7NJyikBIDbQkxtbhzGodSgdohtglKtrC2clAVfnys0WVh3M4tftaaw+lIXZqhLp78GAliEMbBVG10YBuJqkZXeFnCfgLFaVtNOl7EvPZ09aPnvTC9iblk9usXY5eh83E93iAujROIjeTYNoEuItv56EAAm4epZTVM6Kg1ks33+SdUdOUVZpxcPFSJfYBvRsEkSPxoG0CveVrszqOUbAWa0q+aWV5BSXk1NUQXZRBamnS0g5XUJybikpuSWkni6h0qKVaTIoNA31oW2kL20i/WgX5U+bCPt5w7z//vtMmjSJoqIivUsRTuCV67UrXryxtkLnSqrn7e3NpEmTePbZZ/UupVaUVljYkJjN+sRsNh7N5vBJ7TPv7mKgbaQfHaL9aR/tT/sof6IaeMiP8gvZVsDtTs1j6b5MrKoWWharikVVsVpVKiwqJRVmisst2n2FhZJyM3mlleQWV1z0Krv+ni40DPAkOsCThlW3VuG+NA/zsevjUCIiIsjIyNC7DOEk7CngAMLDwy842bIjySosY9OxXHYmnyYhJY+96QVUmK0ABHq50irClxZhPrQI86VluC+NQ7yceTBcjQOuXg6/P5BRwJQ12pH/igJGg4JRUTAYFFyMBrzdjHi6mvByM+Lv4UKEnzv+ni4EerkR4OVKoLcrQd7avyP8PfDzcKmPsuvds88+Ky04IS7C29vbYVpvFxPi487w9hEMbx8BQIXZyuGThexMyWNPah4HMgr54c8kyqtCz2RQaBzsTYtwH5qF+hAX5EWjYC9iA73s+kd+bbPpLkohRB2SfXB2xWyxciKnhAMZBRzMLOBgRiEHMgpIzy+7YLoIP/ezYdcoyIuoBh6E+LoT6utOiI8bLnayy+YybKuLUghhgyTgHEJRuZkT2cUcP+92LLuY46eKKCgz/2X6IG9Xgn3c8fMw4e3mgo+7CW83E97uJrxcjRgMCiaDgkFRtN42g4ICmKt2L5mtKmaL9YK/LVYVs0XFbNUet1Y9fub+/p6xdGzYoLZW2ba6KIUQNsgvWu8KRC3wdjPRJtKPNpF+FzyuqiqnSypJzyslq7CMkwXlnCzQ7rMKyigsM5OeV0pRuVm7lZmpsFivaNkGBUwGA8aqUDQaq+6rdkMZjdp9Xmllba5yjUkLTgghBACVFqs2CPC8gYBmq4qqgotRqQqyc4Fm0Oc4YumiFEII4ZBqHHB2v7dRCCGEuBgJOCGEEA5JAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDkoATQgjhkC55qq7XXnvN9NRTT9VnLUIIIcRlffTRR7FA6quvvvrXE23+j8udizLqo48+qrWihBBCiFpwHGgEnKhuwssFXGrVTOrSmULFlZHtdvVk21092XZXT7bd1bvYtkut0StVVdXtNmnSJFXP5dvrTbabbDvZdvZ1k22nz7bTe5DJazov317Jdrt6su2unmy7qyfb7upd9bar7moCQgghhF3SuwUnhBBC1AkJOCGEEA6pTgNOUZT3FEU5riiKqihKm0tMM0hRlG2KopQrivJeXdZjT2q47f6pKMo+RVESFEXZrijKjfVdpy2q4ba7X1GU3Yqi7FIUZY+iKE/Wd522pibb7bxpmyuKUiKfWU0N33OTFEXJqnrP7VIU5bP6rtMW1fR9pyjK2KrP6t6q+9Dq5l3XLbi5wPVA0mWmOQY8BEyu41rsTU223RYgXlXV9sB44BdFUTzqozgbV5Nt91+gvaqqHYAewLOKorSrj+JsWE22G4qiGIEvq6YXmhptO+AHVVU7VN0er4e67EG1205RlC7AJGCgqqptgF5AfnUzvtxxcNdMVdX1VcVdbprEqmlG1GUt9qaG227peX/uRruUeyA1PUbEQdVw2xWc96cn4AI49Yirmmy3Ki8CCwHvqpvTu4JtJ/5HDbfd08B7qqpmVr2m2nAD2QfnSO4Bjqqq6tThdiUURRmuKMo+tF+Ok1VV3aN3TbauqpV7I/AfvWuxU7dXdY0vUxTlOr2LsSOtgDhFUdYqirJDUZRXlBr8mpCAcwCKovQBXgfu0LsWe6Kq6nxVVVsDzYBxiqI017smW6YoigvwNfCIqqoWveuxQ1OARqqqtkPbJTNPUZRAnWuyFyagHTAQ6AMMAcbV5EXCjlX9CvwJGKGq6iG967FHqqomK4qyBbgJkG14aeFAY2BR1Y9nf0BRFMVXVdW/6VqZHTjTvVb17+WKoqQAbYA1+lVlN5KAX1VVLQfKFUWZB3QFfrjci6QFZ8cURYkHfgHGqKq6Q+967ImiKC3O+3cQ0BeQLsrLUFU1WVXVIFVVY1VVjQU+BL6WcKsZRVEiz/t3ByAW+UFVUzOAQYrGBegPJFT3oro+TOBjRVFSgSjgj6r9HSiKsqhqVAyKovSqmuYZ4GFFUVJluHvNth3wOeABfHne0OO2OpVsM2q47R6uOsRiF7AC+FRV1WU6lWwTarjdxEXUcNu9WTXEPQGtq3fc+a06Z1XDbfczkAXsB3YB+4Bvqp23nKpLCCGEI5IuSiGEEA5JAk4IIYRDkoATQgjhkCTghBBCOCQJOCGEEA5JAk4IIYRDkoATQgjhkCTghBBCOKT/B5HLGe9WQ/tFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_posterior(samples, var_names=['μ'], ref_val=np.log(4));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## To Learn More\n",
    "\n",
    "- Kruschke, J.K. [Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan](https://www.amazon.com/Doing-Bayesian-Data-Analysis-Tutorial/dp/0124058884). 2015. Academic Press / Elsevier. \n",
    "- [Probabilistic Programming and Bayesian Methods for Hackers](http://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/) by Cam Davidson Pilon\n",
    "- [PyMC3 Tutorial Notebooks](https://docs.pymc.io/nb_tutorials/index.html)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
